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Abstract 

Generalized  derivatives  of  fractional  order  are  used 
to  construct  stress -strain  constitutive  relations  for 
viscoelastic  materials,  based  on  the  observed  sinusoidal 
behavior  of  the  materials.  The  non-periodic  behavior  of 
one  material  is  observed  in  the  laboratory  and  compares 
favorably  with  the  non-periodic  behavior  of  the  material 
predicted  by  its  generalized  derivative  constitutive 
relation.  Having  established  that  the  generalized  deriva 
tive  constitutive  relation  is  an  appropriate  mathematical 
model  for  the  general  motion  of  at  least  one  viscoelastic 
material,  the  tools  for  the  analysis  of  structures  of 
engineering  interest  are  put  forward.  In  particular, 
attention  is  focused  on  a finite  element  formulation  of 
and  solutions  to  the  equations  of  motion  for  structures 
containing  elastic  and  viscoelastic  components. 
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APPLICATIONS  OF  GENERALIZED 


DERIVATIVES  TO  VISCOELASTICITY 


I . Introduction 

The  following  investigation  deals  with  constitutive 
relations  employing  generalized  derivatives  that  relate 
stress  and  strain  in  viscoelastic  materials,  and  the  solution 
techniques  for  the  resulting  equations  of  motion  for  structures 
incorporating  viscoelastic  components  to  dampen  vibratory 
motion.  The  use  of  generalized  derivatives  of  fractional 
order  in  stress -strain  constitutive  relations,  first  suggested 
by  Caputo  (Ref  1),  may  be  viewed  as  an  extension  of  the 
standard  model  for  a linear,  viscoelastic  material. 

The  standard  viscoelastic  model  for  a uniaxial  consti- 
tutive relation  is  (Ref  2) 


o(t) 


+ 


K 

l b 


k=l 


k 


F0e(t) 


+ 


d~*  e 
dt^ 


1.1 


The  viscoelastic  constitutive  relation  employing  generalized 
derivatives  of  fractional  order  will  be  taken  to  be 


K 6,  J a . 

o(t)  + l b D k[o(t)]  = E e (t ) + l ED  J[e(t)]  1.2 

k=l  k ° j=l  J 


where  the  generalized  derivative  operator  of  real  order  a is 


defined  by 

5 rTT^rr  ^ TT^a  dT  0 < “ < 1 J-3 

The  generalized  derivative  constitutive  relation,  Eqn  1.2, 
may  be  viewed  as  an  extension  of  the  standard  model,  Eqn  1.1, 
in  the  sense  that  the  derivatives  are  no  longer  limited  to 
being  of  integer  order. 

The  use  of  this  generalized  derivative  constitutive 
relation  in  modeling  the  response  of  viscoelastic  materials 
will  be  seen  to  have  several  advantages  over  present  methods. 

The  major  drawback  of  the  standard  viscoelastic  model, 

Eqn  1.1,  is  that  a large  number  of  terms  are  often  required 
to  describe  a material  adequately.  The  use  of  derivatives 
of  other  than  integer  order  in  the  constitutive  relation  will 
be  seen  to  produce  satisfactory  models  with  very  few  parameters.* 

Because  of  the  large  number  of  terms  required,  the 
standard  model,  Eqn  1.1,  often  becomes  too  cumbersome  to 
manipulate.  Consequently,  an  alternative  known  as  the 
"complex  modulus  method"  has  been  developed.  In  the  complex 
modulus  method,  measured  values  of  E*  (w)  , Eqn  1.4,  are  used 


as  a discrete  approximation  of  the  function  E*(w)  . In  the 

transform  domain,  the  general  viscoelasti*  constitutive 


* In  Section  V,  a generalized  derivative  constitutive  relation 


for  the  elastomer  3M-467  is  presented.  The  relation  charac 
terizes  the  material's  properties  over  four  decades  of 
frequency  with  three  parameters. 


relation  is 


r 


a*  (w)  = E*  (to)  e * (u>)  1.4 

E*(to)  is  measured  for  different  frequencies  of  motion  (Ref  4), 
to^  , which  produces  a set  of  discrete  values  of  the  modulus, 
E*(w^)  , over  the  frequency  range  of  interest.  These  discrete 

values  of  E*(to)  are  substituted  into  the  transformed 
equations  of  motion  of  a viscoelastic  material  to  produce 
values  of  the  transform  of  the  response  at  discrete  frequencies. 
The  inverse  transform  of  the  response  is  evaluated  numerically 
to  produce  the  time  history.  The  major  drawback  of  this  method 
is  the  arduous  task  of  calculating  the  inverse  transform  for 
every  point  in  time  at  which  the  value  of  the  response  is 
required.  The  use  of  the  generalized  derivative  constitutive 
relation  will  be  seen  to  do  away  with  the  need  for  numerical 
approximations  in  the  frequency  domain. 

An  elementary  form  of  the  "complex  modulus"  method, 
obtained  by  representing  the  transform  of  the  modulus  by 


— 


and  often  known  as  "structural  damping"  (Ref  5),  is  valid 
only  for  sinusoidal  stress  and  strain  in  the  material. 

Crandall  has  shown  that  the  response  of  an  harmonic  oscillator 
with  "structural  damping"  for  impulsive  loading  is  non-causal 
(Ref  6);  that  is,  the  time  response  of  the  oscillator  occurs 
before  the  loading. 

Milne  (Ref  7)  has  proposed  several  modifications  of  the 
imaginary  part  of  the  modulus  given  in  Hqn  1.7  to  produce  a 
causal  response.  Ufif ortunately , neither  the  modified  modulus 
nor  the  one  given  in  Eqn  1.7  is  particularly  suitable  for 
transient  (broad-band)  response  of  viscoelastic  materials, 
because  they  do  not  account  for  the  frequency  dependent 
stiffness  typically  encountered.  Caputo  (Ref  8)  observed 
that  a single  term,  generalized  derivative,  constitutive 
relation  of  the  form 

ai 

o(t)  = Exn  1 [ e (t ) ] 0 < < 1 1.7 


i 


i 


I 


produces  frequency -dependent  stiffness  and  damping  and  a 
loss  factor*,  n , that  is  frequency  independent. 


n = tan 


a j tt 


* The  loss  factor  in  a linear  material  is  the  ratio  of 
imaginary  to  the  real  part  of  the  modulus. 


jflRkSf- 


yUKC 


Caputo's  work  with  generalized  derivative  constitutive 
relations  focuses  primarily  on  the  propagation  of  waves  in 
geological  formations.  One  of  Caputo's  earliest  papers 
(Ref  9)  was  on  generalized  independent  loss  factors  or 
equivalently  frequency-independent  resonant  qualities,  Q 
This  work  was  followd  by  a book  (Ref  10)  in  which  Caputo 
dealt  with  the  propagation  of  impulsive  plane  waves  and  the 
free  vibration  of  spherical  strata  using  generalized  deriva- 
tives in  the  constitutive  relations  of  the  media.  In  the  last 
chapter  of  the  book  and  in  later  papers  with  Minardi,  Caputo 
compares  the  generalized  derivative  relations  with  experimental 
observations  of  the  properties  of  some  media:  "some  metals, 
glasses  and  the  earth."  (Refs  11,  12)  In  1974  Caputo 
proposed  a generalized  derivatives  viscoelastic  constitutive 
relation  of  the  form  (Ref  13) 


nDa  n [ e (t ) ] , 0 < a < 1 ; n = 0,1  ,2, . . . 


where  the  generalized  derivative  operator  was  defined  by 


Da+n [x  (t ) ] 


/'  _i dt 

o dtn  ( t - t ) 01 


1.10 


This  definition  differs  slightly  from  the  generalized 
derivative  used  in  this  investigation,  Eqn  1.3.  Caputo 
used  the  relation,  Eqn  1.9,  to  determine  the  response  of  a 
uniformly  driven  infinite  viscoelastic  layer  and  investigated 


the  hysteresis  behavior  of  the  constitutive  relation  (Refs  14 
15).  Recently,  Caputo  suggested,  but  gave  no  application  for 
a constitutive  relation  of  the  form  (Ref  16) 

De[o(t)]  = nDa[ e (t ) ] 1 

which  is  the  harbinger  of  the  general  constitutive  relation, 
Eqn  1.2,  which  is  to  be  used  in  this  study. 

All  of  Caputo 's  work  rests  on  a continuum  formulation 
of  the  equations  of  motion.  Unfortunately,  a continuum 
formulation  of  the  equations  of  motion  for  many  structures 
of  engineering  interest  is  not  practical.  As  a result,  a 
discrete  formulation  of  the  equations  of  motion  is  adopted, 
based  on  assumed  displacement,  finite-element  methods.  A 
solution  technique  developed  for  the  motion  of  structures 
incorporating  viscoelastic  materials  modeled  with  generalized 
derivatives  is  developed  as  an  extension  of  the  solution 
technique  developed  by  Foss  for  non-proportional  viscous 
damping  (Ref  17). 

In  the  sections  to  follow,  constitutive  relations  using 
generalized  derivatives  are  developed,  and  their  applications 
and  limitations  are  considered.  Some  experimental  results 
are  presented  and  found  to  show  that  a constitutive  relation, 
with  parameters  determined  from  the  response  to  sinusoidal 
loading,  predicts  very  well  the  response  of  one  typical 


elastomeric  material  to  an  impact  loading.  Finally,  the 
formulation  of  multi-degree  of  freedom  systems,  necessary 
for  large-scale  structural  analysis,  is  considered.  Special 
solution  methods,  necessary  for  such  applications,  are 
developed. 
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A Brief  Overview  of  Generalized  Derivatives 


Before  proceeding  to  construct  generalized  derivative 
constitutive  relations,  it  is  appropriate  to  introduce  the 
properties  of  generalized  derivatives  relevant  to  the  follow- 
ing investigations.  Of  particular  interest  are  the  form  of 
the  Laplace  and  Fourier  transforms  of  generalized  derivatives, 
and  the  results  of  repeated  differentiation  of  fractional 
order . 

First  and  foremost,  the  generalized  derivative  is  a 
linear  operator. 

Da[x1  (t ) + x2(t)]  = Da[x1 (t) ] + Da[x2(t)]  2.1 


5 


This  property  follows  directly  from  the  definition,  Eq  1.3. 


Da[x(t)] 


r (1  -a)  cTT 


d ft  xCt) 


dt 


o ( t - T ) 


1.3 


To  put  the  definition  into  a form  in  which  the  calculation 
of  its  Laplace  transform  is  straightforward,  one  first  performs 
a change  of  variable 


t = t - n 


which  results  in 


2.2 


D“[x(t)]  - iTTnrr  A l d" 


2.3 


Using  Leibnitz's  rule  to  differentiate  the  integral  produces 


D°[x(t)]  ' rnrrr  / A A dn  * rrr'°rt»  2-4 

After  taking  the  Laplace  transform  of  Eq  2.3,  the  transform 
of  the  generalized  derivative  of  order  a of  the  function 
x(t)  is  seen  to  be 


L[Da[x(t)]] 


1 - a 


(s  L[x  (t ) ] - x (0 ) ) 


+ x(0) 


1 - a 


2.5 


which  simplifies  to 


L[Da[x(t)]]  = s aL[x (t ) ] 


2.6 


where 


L[x (t ) ] 


/ x(t)c~st  dt 


2.7 


Notice  that  the  Laplace  transform  of  a generalized  derivative 
of  order  a of  a function  is  equal  to  sa  times  the  transform 
of  the  function. 

Under  certain  conditions  a similar  property  of  generalized 
derivatives  is  true  for  Fourier  transforms. 


F[Da[x(t)]]  = (iu»)a  F[x  (t ) ] 


2.8 


- 


— 


4 


where 


oo 

F [x  (t ) ] = j x (t ) e” lwt  dt 

- 00 

The  conditions  are, first,  that 


2.9 


x(t)  = 0 for  t < 0 


2.10 


in  which  case  the  Fourier  transform  becomes 

00  . 

F[x  (t ) ] = / x(t)e'lwt  dt  2.11 

o 

and,  second,  that  the  integral  in  Eq  2.11  exists.  Note  the 
parallel  form  of  Eqns  2.6  and  2.8.  Both  relations  were  used 
by  Caputo  (Ref  18). 

A useful  property  of  generalized  derivatives  is  that 
the  generalized  derivative  of  order  , of  the  generalized 

derivatives  of  order  of  a function  is  the  generalized 

derivative  of  order  + 0^  of  the  function.  Again  using 

operator  notation,  the  property  is 


[D“2[x(t)]] 


+ a . 


[x(t)  ] 


2.12 


Notice  that  the  definition  of  generalized  differentiation 
given  in  Eq  1.3  is  restricted  to  fractional  order  a less 
than  one.  If  a is  one  or  greater  in  Eq  1.3,  the  integral 


10 


m'  ."/V 


The  definition  of  a 


* 

contains  a non- integrable  singularity  . 
generalized  derivative  of  order  B , where  B > 1 and 
B = m+a  where  m is  the  largest  integer  not  exceeding  B , 
is 


Dm+a  [x(t)] 


X (t) 
(t-T)0 


dt 


2.13 


* On  the  other  hand,  if  a is  zero  in  Eq  1.3,  the  relation  is 
clearly  valid  and  follows  from  the  fundamental  theorem  of 
calculus . 


r ( 

H 


li 


In  this  chapter  the  basic  generalized  derivative, 
viscoelastic  constitutive  relation  is  presented  and  some 
aspects  of  its  behavior  are  established.  The  basic  general- 
ized derivative  constitutive  relation  is 


o(t)  = XDa[e(t)]  0 < a < 1 


3. 


or 


a (t) 


t 

I 


o 


e (n) 
(t-n)a 


dn 


3. 


Since  the  generalized  derivative  is  a linear  operator, 
this  relation  is  suitable  only  for  the  linear  approximation 
of  a material's  properties.  This  linear  constitutive  relation 
satisfies  many  of  the  presently  accepted  constraints  on 
viscoelastic  constitutive  relations. 

In  particular,  it  represents  a material  with  fading 
memory.  To  demonstrate  this  claim,  it  is  necessary  to  put 
the  constitutive  relation  in  a different  form  using  a change 
in  variable 


T|  = t - T 


3. 


12 


and  again  using  Leibnitz's  rule  to  differentiate  the  resulting 
integral  produces 


o(t) 


rl  1 3 e(t-t)  , 

TT  [ ~ at  dT  + 


Xe(0) 

r(l-a)t' 


(t)  = / G(T)e(t-T)dT  + G(t)e(0) 


where 


G(t) 


r(l-a)ta 


The  constitutive  relation,  Eq  3.5,  represents  a visco- 
elastic material  with  a fading  memory.  A material  is  said 
to  have  a fading  memory  if  its  relaxation  modulus,  G(t)  , 
goes  to  zero  monotonical ly  as  t increases  (Ref  20).  Notice 
that  G(t)  in  Eq  3.6  does  in  fact  go  to  zero  monotonically . 

Since  the  material  has  a memory,  the  value  of  the  stress 
at  time  t is  dependent  on  the  entire  strain  history  up 
until  time  t . To  insure  that  the  constitutive  relation 
produces  a stress  that  is  dependent  on  the  entire  strain 
history  up  until  t , time  zero  must  be  chosen  before  or 
at  the  onset  of  the  initial  strain. 

Consequent ly , 


0 for  t < 0 
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and  the  only  way  that 


e(0)  t 0 


is  if  the  strain  history  is  discontinuous  at  t = 0 
According  to  Gurtin  and  Sternberg  (Ref  21),  the  constitutive 
relation  as  shown  in  Eq  3.5  is  in  the  correct  form  to  handle 
discontinuous  strain  histories.  Notice  that  a step  discontin- 
uity in  the  strain  history  at  t = 0 produces  a stress 
* 

history  that  is  singular  at  t = 0 

When  the  strain  history  is  a continuous  function  of 
time,  and  zero  for  negative  times,  the  constitutive  relation 
given  in  Eq  3.5  reduces  to 


(t)  = / G(T)e(t-T)dT  , G(t) 


r(l-a)t 


This  relation  satisfies  three  out  of  four  of  Pipkin's  restrictions 
on  viscoelastic  constitutive  relations  (Ref  22).  The  first 
restriction,  that  the  stress  be  an  odd  functional  of  strain 
rate,  is  satisfied.  The  second  restriction,  that  G(t)  go  to 
zero  as  time  increases,  is  also  satisfied,  which  is  in  keeping 
with  the  fading  memory  property.  The  third  restriction,  also 
satisfied,  is  that  the  kernel,  the  relaxation  modulus  G(t)  , 
be  a function  and  not  a distribution.  Pipkin's  fourth 


The  Voigt  viscoelastic  model  displays  this  same  property 


m$K> 


restriction,  not  satisfied  by  the  generalized  derivative 
constitutive  relation,  is  that  G(t)  be  of  negative 
exponential  order. 


IV . Generalized  Derivative  Constitutive 
Relations  for  Viscoelastic  Materials 


The  task  at  hand  is  to  use  the  basic  generalized 
derivative  constitutive  relation,  just  presented,  as  the 
building  block  for  constitutive  relations  that  model  the 
frequency  dependent  moduli  of  viscoelastic  materials.  The 
moduli  of  viscoelastic  materials  are  complex  numbers  where 
the  real  and  imaginary  parts  are  functions  of  the  frequency 
of  motion. 


A*  (oo) 

= X " (10)  + 

i X ' ' (u>) 

4.1 

P*  (w) 

= p ' (w)  + 

ip  " (w) 

4.2 
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at  frequency  u>Q  relate  sinusoidal  stress  and  strain  of 
frequency  in  the  material. 


o ft) 
mn v 


6 x * (u  )e  exp[iw  t]  + 2y * (u  )e  „ exp[iw  t]  4.S 
mn  o o 1 o ^ o mn  r o 
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Consequently,  one  can  measure  the  values  of  y*(u>)  and 
X*  (w)  at  discrete  frequencies  of  sinusoidal  motion.  As  a 
result,  the  frequency  dependence  of  the  moduli  can  be 
determined  experimentally. 

Typically,  a viscoelastic  material  at  constant,  uniform 
temperature  has  moduli  that  vary  with  the  frequency  of  motion 
as  indicated  in  Pig  4.1  (Ref  23).  At  low  frequencies  (the 
rubbery  region)  the  real  part  of  the  modulus  is  relatively 
constant,  while  the  imaginary  part  of  the  modulus  increases 
with  increasing  frequency.  At  intermediate  frequencies  (the 
transition  region)  both  the  real  and  imaginary  parts  of  the 
modulus  increase  with  increasing  frequency  and  the  rate  of 
increase  of  the  real  part  slowly  overtakes  the  rate  of 
increase  of  the  imaginary  part.  At  high  frequencies  (the 
glassy  region)  the  imaginary  part  of  the  modulus  decreases 
with  increasing  frequency,  and  the  real  part  of  the  modulus 
is  relatively  constant. 

The  generalized  derivative  constitutive  relations 
presented  here  are  of  two  types.  The  first  type  are  those 
relations  intended  to  model  the  viscoelastic  behavior  of 
the  material  in  the  rubbery  and  transition  regions.  For 
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brevity,  this  type  of  model  is  referred  to  as  the  RT  model. 
The  second  type  are  those  relations  intended  to  model  the 
behavior  of  the  material  in  the  rubbery,  transition  and 
glassy  regions.  This  type  of  relation  is  referred  to  as 
the  RTG  model.  Since  the  RTG  model  may  be  viewed  as  a 
generalization  of  the  RT  model,  the  RT  model  is  considered 
first . 

The  RT  model  for  an  isotropic,  homogeneous,  linear 
viscoelastic  material  is 
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To  establish  the  frequency  dependence  of  the  moduli 
in  the  RT  model,  one  takes  the  Fourier  transform  of  the 


and 


I v*  ^ 

i = l a 


4.14 


the  moduli  in  the  RT  model,  Eqns  4.11  and  4.12,  have  essentially 
a constant  real  part  and  an  imaginary  part  that  increases  with 
increasing  frequency,  similar  to  the  properties  of  a visco- 
elastic material  observed  in  the  rubbery  region,  Fig  4.1. 

At  intermediate  frequencies  of  motion,  defined  by 


a 

Y X . ui 

= 1 J 
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4.15 


and 

A 

1 V - 1 

— - 1 • 

O 1=1 

the  real  and  imaginary  parts  of  the  moduli  n the  RT  model 
are  increasing  with  frequency,  similar  to  the  properties  of 
a viscoelastic  material  observed  in  the  transition  region. 

It  is  evident  that  the  RT  model  does  not  properly 
account  for  the  properties  of  a viscoelastic  material  at 
high  frequencies,  defined  by 
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The  real  and  imaginary  parts  of  the  modulus  are  predicted  by 
the  RT  model  to  increase  indefinitely  with  frequency.  In  the 
glassy  region,  however,  the  real  part  of  the  modulus  of 
viscoelastic  materials  is  typically  constant,  Fig  4.1,  and 
the  imaginary  part  is  decreasing  with  increasing  frequency. 
This  discrepancy  motivates  the  construction  of  a more  complex 
model  which  accounts  for  material  properties  in  the  glassy 
region,  the  RTG  model. 

The  RTG  model  for  an  isotropic,  homogeneous,  linear, 
viscoelastic  material  is  defined  to  be 
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Again,  the  frequency  dependence  of  the  moduli  in  the  model  is 
observed  by  taking  the  Fourier  transform  of  the  constitutive 
relation.  After  some  algebraic  manipulation  of  the  transformed 
relation,  the  result  is 
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where 


A (w) 


J a. 

(a  + l \.(i«o)  D) 
° 3 = 1 J 

( 1 * l b (ito)  k) 

p=l  P 


4.21 


U*  (w) 


(w0  + l *) 

0 t=l  1 

p r~ 

Cl  * l b (ia>)  P) 

p=l  p 


4.22 


0 < a . , a „ , B,,  B < 1 
3 l k p 


4.23 


and 


Xo’Xj  ’yo,lVak’bp  > 0 


4 . 24 


2 3 


hA  > ••  * >* 


K> 

4 


If  the  parameters 
that 


a,  and  b are  chosen  to  be  small,  such 
k p ’ 


K 

J 

a . 

I a.  (iio)  k 
k=l  k 

< < 

X + l X . 

° j-1  J 

CiuO  J 

and 


P 

L a . 

£ b (iw)  P 

i P 

p=l  ** 

< < 

v*o  + l W,(iw) 

0 *=1  * 

4.26 


for  frequencies  in  the  rubbery  and  transition  regions  of 
the  material,  the  RTG  model  behaves  like  the  RT  model  in  the 
ruhbery  and  transition  regions  of  the  material. 

It  is,  however,  the  presence  of  the  terms  and  b^ 

which  enable  the  RTG  model  to  account  for  the  properties  of 
viscoelastic  materials  at  high  frequencies;  i.e.,  in  the 
glassy  region.  Note  that,  if  the  largest  values  of  and 

* A 

gj,  are  the  same,  and  if  the  largest  values  of  and  g 

* * 

are  the  same,  then  at  high  frequencies  X (w)  and  u (w) 
have  real  parts  that  become  constant  and  imaginary  parts  that 
decrease  with  increasing  frequencies , as  is  characteristic  of 
a viscoelastic  material  in  the  glassy  region. 

An  important  property  of  the  RT  and  RTG  models  is  that 
they  satisfy  the  "elastic-viscoelastic  correspondence  princi- 
ple." (Ref  4)  The  correspondence  principle  states  that  the 
Laplace  transform  of  the  stress  response  of  a viscoelastic 
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material  can  be  constructed  from  the  Laplace  transform  of 
the  response  of  an  elastic  material  by  replacing  the  elastic 
constants,  \ and  p , in  the  elastic  response  by  the 
Laplace  transforms  of  viscoelastic  moduli,  A*(s)  and 
(s)  . The  principle  holds  when  the  transform  of  the 

elastic  stress -strain  constitutive  relation 
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can  be  used  to  construct  the  transform  of  the  viscoelastic 
stress  - strain  constitutive  relation  by  replacing  the  elastic 
constants  with  the  transforms  of  the  viscoelastic  moduli. 
Thus,  the  viscoelastic  constitutive  relation  must  be  of  the 
form 
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The  Laplace  transforms  of  the  RT  and  the  RTG  models  are 
of  the  general  form  given  in  Eq  4.28.  The  Laplace  transform 
of  the  RTG  model  is 
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and  the  Laplace  transform  of  the  RT  model  is 


4.30 


Another  important  property  of  the  RT  and  RTG  models  is 

that  they  can  be  constructed  to  be  causal  in  the  sense  that 

the  response  (stress)  does  not  occur  before  the  input  (strain). 

The  stress  response  is  zero  for  negative  time  if  its  Laplace 

transform  is  analytic  in  the  right  half  s plane.  This  condition 

on  the  transform  of  the  stress  is  met,  for  the  class  of  strain 

histories  having  transforms  that  are  analytic  in  the  right 

a-  a» 

half  s plane,  when  the  branch  cuts  of  s J , s , s 
BD 

and  s F are  along  the  negative,  real  s axis,  and  y*(s) 
and  A*(s)  have  no  zeros  in  the  right  half  s plane.  Since 
the  stress  is  zero  for  negative  time,  the  stress  cannot 
anticipate  the  strain  that  begins  at  time  zero. 

Laplace  transforms  are  also  useful  in  determining  the 
hysteresis  predicted  by  the  RT  and  RTG  models.  For  a 
sinusoidal  strain  history  of  frequency  starting  at  time 

zero,  the  transform  of  the  stress  history  is 
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Evaluating  the  inverse  transform  in  the  same  manner  as  in 
Section  X,  the  stress  time  history  is  found  to  be 
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As  time  increases,  the  integral  term  in  Eq  4.33  goes  to  zero 
and  the  sinusoidal  terms  dominate  the  stress  time  history. 
So,  for  a sinusoidal  strain  history,  the  stress  eventually 
becomes  sinusoidal  as  well. 

Using  Euler's  formula,  numerous  trigonometric  manipu- 
lations, and  the  observation  that  the  two  sinusoidal  terms 
are  conjugates,  the  stress  for  large  time  may  be  evaluated 
as  components  in-phase  with  the  strain,  proportional  to 
sinwot  , and  components  out -of -phase , proportional  to 
cosioQt  . The  resulting  expression  for  the  stress  using  the 
RTG  model  is 
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and 
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A more  compact  form  of  the  expression  for  the  stress  is 
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Under  conditions  of  uniaxial  stress  and  strain,  Eqns  4.31 
and  4.37  may  be  recognized  as  the  parametric  equations  of 
an  ellipse;  thus,  it  is  clear  that  the  RTG  model  predicts 
the  existence  of  a hysteresis  loop  and  the  loop  is  elliptical. 

The  loss  factor  associated  with  the  hysteresis  loop,  the 
ratio  of  the  energy  dissipated  during  a cycle,  D,  to  the  peak 
strain  energy  stored  during  a cycle,  U , is  a parameter 
often  used  to  characterize  the  ability  of  viscoelastic 
materials  to  dampen  vibratory  motion. 
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The  energy  dissipated  per  cycle  is 
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where 


— is  the  period  of  the  motion.  Using  Eqns  4.31  and 
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4.37  to  evaluate  dissipated  energy  yields 
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The  peak  energy  stored  during  a cycle  is 
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where  o are  the  stresses  in-phase  with  their  respective 

strains,  t sin«  t . Again  using  Eqns  4.32  and  4.33,  the 
mn  o 

o 

peak  energy  stored  during  a cycle  is 
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The  resulting  loss  factor  is  seen  to  be 
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Notice  that  the  form  of  the  expressions  for  the  hysteresis 
behavior  and  loss  factor  of  the  RT  model  is  identical  to  that 
of  the  RTG  model.  This  follows  from  the  observation  that  the 
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RT  model  is  a special  case  of  the  RTG  model  for  = 

k = 1 , 2 , . . . , K , and  b^  = 0 , p = 1 , 2 , . . . , P 
In  summary,  RT  and  RTG  models  satisfy  the  elastic 
viscoelastic  correspondence  principle.  Conditions  necessary 
to  insure  that  the  stress  does  not  anticipate  the  strain  have 
been  developed.  In  addition,  both  models  predict  the  existence 
of  stress -strain  hysteresis  effects  and  the  resulting  hysteresis 
loops  are  elliptical.  Most  significantly,  the  models  predict 
moduli  which  have  the  same  frequency  dependence  as  is  observed 
in  frequency -dependent  moduli  in  typical  viscoelastic  materials.* 
The  outstanding  question  is  whether  or  not  the  parameters 
of  the  RT  and  RTG  models  can  be  chosen  to  describe  accurately 
the  properties  of  a particular  viscoelastic  material.  The 
construction  of  the  RT  model  for  the  elastomer  3M-467  is 
the  topic  of  the  following  section. 
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* In  addition,  a generalized  derivative  constitutive  relation  occurs 
in  Newtonian,  viscous  fluids  as  demonstrated  in  Appendix  B. 


V . The  RT  Model  for  the  Elastomer  5M-467* 

The  first  material  examined  for  the  possible  application 
of  generalized  derivative  constitutive  relations  was  the 
adhesive  tape  3M-467.  The  tape  was  chosen  as  a prime  candi- 
date because  of  its  viscoelastic  mechanical  properties,  its 
linear  response  in  shear  for  engineering  shear  strains  up 
to  1,  its  growing  applications  in  mechanical  damping,  and 
the  fact  that  sufficient  data  on  its  mechanical  properties 
were  available. 

The  proposed  uniaxial  shear  RT  model  for  3M-467  is 

V/1’  ‘ 2 <"o  * "lDai>  Sm(t)  ' m * n 5 


The  parameters  of  the  model,  ’ yl  ’ an<^  al  » 

are  chosen  such  that  the  sinusoidal,  steady-state  response 
of  the  model  closely  approximates  the  sinusoidal,  steady- 
state  response  of  the  material  observed  experimentally.  For 
sinusoidal  strain 
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the  RT  model  generates  stresses  of  the  form 
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for  t large  enough  for  the  transients  to  have  died  out. 
The  frequency  dependent  shear  modulus,  u(w)  , is  seen 
to  be 
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Figure  5.1  displays  the  good  agreement  between  the  experimen- 
tally observed  mechanical  properties  of  3M-467  at  75°F  and 
the  RT  model  using  the  values  of  the  parameters  given  above.* 
The  parameters  of  the  RT  model  were  determined  in  an 
iterative  manner.  Initial  guesses  of  the  parameters  were 
made,  and  the  resulting  frequency  dependent  shear  modulus 


* The  mechanical  properties  of  3M-467  were  provided  by  the 
U.S.  Air  Force  Materials  Laboratory,  Wright  - Patterson 
Air  Force  Base,  Ohio. 
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Figure  5.1.  A Comparison  of  the  Material  Properties  of  34-467 
and  the  Properties  Predicted  by  its  RT  Model 


was  compared  to  the  observed  modulus.  Successive  guesses 
of  the  parameters  were  made  to  match  the  slopes  and 
asymptotes  of  the  model  to  those  of  the  observed  properties 
until  an  acceptable  fit  was  obtained. 

Although  the  parameters  of  the  RT  model  are  based  on 
the  sinusoidal  response  of  the  material  at  75°F,  the  model 
can  be  used  for  non-periodic  strain  histories.  To  demonstrat 
the  ability  of  the  model  to  portray  accurately  the  behavior 
of  the  material  when  undergoing  non-periodic  motion,  the 
response  of  the  material  as  predicted  by  the  RT  model  is 
compared-  to  the  experimentally  observed  response  of  the 
material  at  75°F. 

In  particular,  the  behavior  of  3M-467  was  observed  when 
the  material  was  used  as  a viscoelastic  spring  in  a simple 
oscillator  undergoing  non-periodic  motion.*  The  visco- 
elastic spring  was  two  pads  of  3M-467  that  underwent'  shear 
strain  during  the  motion  of  the  oscillator.  Each  pad  of 
3M-467  was  made  by  laminating  2 mil  layers  of  3M-467.  Air 
entrapped  between  the  layers  of  the  pad  was  removed  by 
pressing  the  layers  together  with  a 5 lb.  weight  for  48  hours 

The  two  other  components  of  the  oscillator  are  the  mass 
and  the  support  structure,  Fig  5.2.  The  mass  for  the 
oscillator  is  a metal  cube  sandwiched  between  the  two  visco- 
elastic pads,  Fig  5.3.  Each  pad  is  attached  to  an 


* A schematic  of  the  viscoelastic  spring  appears  in  Fig  5.3. 


Figure  5.2.  A Schematic  of  the  Oscillator  and  its  Support  Structure 


Figure  5.3.  Expanded  Schematic  of  the  Oscillator  and  Support 
Structure  Not  Displaying  the  Steel  Foundation 


aluminum  brace.  Both  braces  are  glued  to  an  aluminum  base 
which,  in  turn,  is  glued  to  a steel  foundation  as  shown 
in  Fig  5.2. 

The  specific  objective  of  the  experiment  was  to 
determine  the  acceleration  transfer  function  of  the  oscilla- 
tor. The  acceleration  transfer  function  is  the  ratio  of 
the  transform  of  the  acceleration  time  history  to  the  trans- 
form of  the  input  force  history.  The  force  time  history, 
measured  by  a Wilcoxon  Z-ll  impedance  head,  was  sampled 
at  2 x 104  measurements  per  second  and  all  frequencies  above 
8 x 103  Hz  are  filtered  out.  The  mass  of  the  oscillator  was 
tapped  with  a Wilcoxon  Z-ll  impedance  head  to  produce  impul- 
sive loading.  The  force  time  history  measured  by  the  impedance 
head  and  the  resulting  acceleration  time  history,  measured  by 
an  Endevco  accelerometer.  Model  2217,  were  also  sampled  at 
2 x 104  measurements  per  second  where,  as  before,  all 
frequencies  above  8 x 103  Hz  were  filtered  out.  The  trans- 
forms of  the  time  histories  were  calculated  using  the  "fast 
Fourier  transform”  routines  of  the  Hewlett  Packard  System  5451B. 

This  experimentally  determined  transfer  function  is 
compared  to  the  analytically  predicted  transfer  function 
based  on  the  equations  of  motion  of  the  oscillator  and  the 
RT  model  for  the  viscoelastic  pads.  The  force-displacement 
relation  for  the  two  pads  based  on  the  RT  model,  Eqn  5.1,  is 


5.8 


fP(t) 


) x(t) 


where  A is  an  area  of  contact  with  the  mass  for  each  pad, 

6 is  the  thickness  of  the  pad,  f (t)  is  the  total  force 
acting  on  the  faces  of  the  pads  in  contact  with  the  mass, 
and  x(t)  is  the  displacement  of  the  face  of  the  pad  in 
contact  with  the  mass.  This  force  displacement  relation  is 
based  on  the  assumption  that  the  displacement  in  the  pad 
varies  linearly  between  the  support  wall  and  the  mass  of 
the  oscillator.* 

The  resulting  equation  of  motion  for  the  oscillator  is 

A 

f(t)  = mx  (t)  + ^ (uo  + ^D*1)  x(t)  5.9 


Taking  the  Fourier  transform  of  the  equations  of  motion  and 
determining  the  acceleration  transfer  function  produces 
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5.10 


A comparison  of  the  experimentally  determined  and 
analytically  predicted  transfer  functions  for  five  oscillators 
with  various  masses  and  viscoelastic  spring  stiffnesses  is 
presented  in  Figs  5.4  through  5.13.**  Each  transfer 


* A finite  element  analysis  of  the  viscoelastic  pad  verifies 
this  assumption  to  be  valid  for  the  frequency  range  of  the 
tests,  0 to  5 x 103  Hz. 

**  The  relevant  parameters  for  each  of  the  five  oscillators 
arc  given  in  Table  5.1. 
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Includes  the  mass  of  the  accelerometer. 


Frequency  (H 


Frequency  (Hz) 


Frequency  (He) 


qucr.c 


Frequency  (He) 


Frequency  (Hz) 


rcqucncy  (Hz) 


function  is  displayed  in  terms  of  its  magnitude  and  phase. 

The  agreement  between  the  observed  and  predicted  transfer 
functions  is  very  good. 

For  comparison,  the  calculated  transfer  function 

* 

based  on  a Voigt  viscoelastic  model  of  3M-467 

0mn(t)  * 2 (lVnm(t)  * * m * n 

is  also  given  in  Figs  5.4  through  5.13.  The  parameters  of 
the  Voigt  model,  and  , are  chosen  to  match  the 

properties  of  3M-467  at  103  Hz. 

Pq  = 630  lb/in2 

= .113  lb-sec/in2 

Note  that  for  some  of  the  oscillators,  the  Voigt  model 
and  the  RT  model  both  generate  transfer  functions  that  agree 
reasonably  well  with  the  observed  transfer  functions.  However, 
for  those  oscillators  having  the  peak  magnitude  of  acceleration 
response  at  higher  frequencies.  Fig  5.12  for  example,  the 
transfer  function  based  on  the  RT  model  is  clearly  in  better 
agreement  with  the  measured  transfer  function  than  the 
transfer  function  based  on  the  Voigt  model.  In  addition, 
the  phase  of  the  observed  transfer  functions  is  consistently 
modeled  more  accurately  by  the  phase  of  the  transfer  functions 

\ 

\ 


5.12 

5.13 
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Figure  5.14.  A Comparison  of  the  Material  Response  of  3M-467  and  the 
Properties  Predicted  by  a Voigt  Model. 


calculated  using  the  RT  model.  These  results  follow  directly 
from  the  fact  that  the  RT  model  accounts  for  the  observed 
properties  of  3M-467  over  the  entire  frequency  range  of 
interest,  102  Hz  to  5 x 103,  whereas  the  Voigt  model  accounts 
for  the  observed  properties  of  3M-467  only  in  the  neighborhood 
of  103  Hz.  This  is  clearly  seen  by  comparing  Figs  5.1  and  5.14. 

If  one  attempts  to  duplicate  the  results  presented  here, 
one  should  be  aware  that  the  mechanical  properties  of  3M-467 
are  strongly  dependent  on  the  water  present  in  the  material . 
Figure  5.15  shows  the  variation  of  the  real  part  of  the 
modulus  with  relative  humidity.  Changes  in  the  imaginary 
part  of  the  modulus  with  relative  humidity  are  roughly 
proportional  to  changes  in  the  real  part.  Hence,  the  loss 
factor,  the  ratio  of  imaginary  part  to  real  part  of  the 
modulus,  is  relatively  insensitive  to  changes  in  relative 
humidity,  as  seen  in  Fig  5.16. 

The  pads  of  3M-467  used  in  this  experiment  were  fabri- 
cated under  conditions  of  40%  relative  humidity  at  room 
temperature.  However,  the  pads  were  kept  covered  during  the 
time  between  fabrication  and  installation  into  the  test  setup. 

In  summary,  the  RT  model  for  3M-467  is  capable  of 
accurately  predicting  the  non-periodic  response  of  the  material 
over  several  decades  of  frequency,  and  is  superior  to  a Voigt 
model  of  the  material.  Consequently,  the  RT  model,  having 
parameters  based  on  the  sinusoidal  steady  motion  of  the  mater- 
ial at  numerous  frequencies,  is  capable  of  predicting  the 
response  of  the  material  to  impulse-like,  short  duration 
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Figure  5.15.  The  Variations  of  Young's  Modulus 

of  3M-467  with  Chano.es  in  Relative 
Humidity. 


Figure  5.16.  The  Variations  of  Loss  Factor  of  3M-467 
with  Changes  in  Relative  Humidity 


loading.  Therefore,  one  can  conclude  that  the  RT  model 
can  accurately  predict  the  general  response  of  the  material 
within  the  frequency  range  of  the  model. 


VI . A Finite  Element  Formulation 
of  the  Equations  of  Motion 


Having  established  that  a very  elementary  formulation 
of  the  equation  of  motion  for  an  elastomer-damped  oscillator 
produced  excellent  agreement  with  experimental  observation,  it 
is  appropriate  at  this  point  to  put  forward  the  tools  required 
in  the  analysis  of  more  complicated  structures  of  engineering 
interest.  In  particular,  the  development  focuses  on  the 
analysis-  of  structures  having  both  elastic  and  viscoelastic 
component  s . 

A continuum  formulation  of  the  equations  of  motion  for 
such  structures  is  impractical  because  of  the  resulting 
complexity  of  the  formulation  for  most  structures  with  complex 
geometry  and  varying  material  properties.  As  a result,  a 
finite  element  formulation  of  the  equations  of  motion  is 
adopted . 

The  cornerstone  of  the  finite  element  approach  is  the 
construction  of  the  stiffness  matrices  for  each  of  the 
finite  elements  in  the  structure.  The  stiffness  matrices  for 
the  elastic  finite  elements  of  structure  are  constructed  in 
the  normal  fashion  using  assumed  displacement  methods  or 
assumed  stress  methods,  etc. 

The  formulation  of  the  stiffness  matrices  for  the  finite 
elements  in  the  viscoelastic  components,  however,  is  limited 
to  those  methods  that  do  not  constrain  the  stresses  in  each 
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finite  element  to  be  in  equilibrium  with  the  forces  at  the 
nodes  of  the  element.  The  assumed  stress  method,  in  parti- 
cular, is  based  on  this  constraint  (Ref  24).  As  a result, 
the  time  dependence  of  the  stresses  in  the  element  is  predi- 
cated on  time  dependence  of  the  nodal  forces.  However,  this 
contradicts  the  fundamental  nature  of  the  generalized 
derivative  models  in  which  the  time  dependence  of  the  stresses 
is  predicated  on  the  time  dependence  of  the  strain  histories. 

Consequently,  the  assumed  displacement  method  is  adopted 
to  formulate  the  stiffness  matrix  of  the  viscoelastic  finite 
element  (Ref  25).  In  an  assumed  displacement  element,  the 
displacements  within  the  element  are  assumed  functions  of 
the  nodal  displacements. 

The  stiffness  matrix  of  the  viscoelastic  finite  element 
is  constructed  using  the  elastic-viscoelastic  correspondence 
principle.  The  stiffness  matrix  is  first  formulated  as 
though  the  material  were  elastic.  The  stiffness  matrix  is 
then  separated  into  two  matrices,  one  matrix  containing 
those  elements  proportional  to  the  elastic  constant  X , 
and  the  other  matrix  containing  those  elements  proportional 
to  the  elastic  constant  u 

[K  ] =A[K]+n[K]  6.1 

eee 

At  this  point  the  transforms  of  the  moduli,  p*(s)  and 
X*(s)  , from  either  the  RT  or  the  RTG  models,  are  substi- 

tuted in  place  of  the  elastic  constants,  v and  X 
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The  result  is  the  viscostiffness  matrix  of  the  finite 
element,  [Ke(s)] 

[Ke(s)]  = X*(s)[r]  + y*(s)tK'-] 


Specifically,  the  viscostiffness  matrix  for  a finite  element 
in  which  the  RTG  constitutive  relation  is  used  to  model  the 
material  is 


[Ms)]  - 


J a . 
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6.3 


The  viscostiffness  matrix  of  the  finite  element  relates 
the  nodal  forces,  {F(s)}  , and  nodal  displacements,  (X(s)}  , 

as  shown  below 


{F  (s)  } = [Kg  (s  ) ] {X(s)> 


and  the  viscostiffness  matrix  of  a viscoelastic  structural 
component  is  constructed  from  the  viscostiffness  matrices  of 
the  elements  within  the  component  in  the  normal  manner.  The 
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viscostiffness  matrices  of  the  R viscoelastic  structural 
components 

{F(s)}r  = [K(s)]r  {X(s)}r  r = 1,2,3 R 

and  the  stiffness  matrices  of  the  Q elastic  structural 
components 

{F(s)}q  = [K]q  (X (s ) }q  q = 1,2,3,.  . . ,Q 

are  used  to  construct  the  stiffness  matrix  of  the  total 
structure,  again  in  the  normal  manner. 

The  stiffness  matrix  of  the  total  structures,  [K(s)J 
and  the  mass  matrix  of  the  total  structure  are  now  used  to 
construct  the  Laplace  transform  of  the  equations  of  motion 
of  the  structure 


s2[M]  (X  (s ) } + [ K (s ) ] { X ( s j } = {F  (s ) } 

Since  some  of  the  elements  in  [ K ( s ) ] are  functions  of  s 
decoupling  the  equations  of  motion,  Eqn  6.7,  to  obtain 
solutions  is  more  complicated  than  decoupling  the  equations 
of  motion  of  a completely  elastic  structure  where  the 
stiffness  matrix  has  constant  elements.  Finding  solutions 
to  Eqn  6.7  is  the  topic  of  the  next  section. 


VII.  The  Solution  of  the  Discrete 


Equations  of  Motion 

The  task  at  hand  is  the  solution  of  the  equations  of 
motion  which  resulted  from  the  finite  element  formulation, 

Eqn  6.7.  A form  of  modal  analysis  is  adopted  where  the  mode 
shapes,  eigenvectors,  of  the  equations  of  motion  are  used 
to  construct  an  orthogonal  transformation  of  the  variables 
that  decouples  the  equations  of  motion.  The  decoupled 
equations  of  motion  are  then  used  to  determine  the  compon- 
ents of  the  structure's  response  and  a general  form  of  the 
solution  to  the  equations  of  motion  is  derived. 

Throughout  this  development,  the  viscoelastic 
components  of  the  structure  are  described  by  their  respec- 
tive RTG  models.  Since  the  RT  model  is  a simplified  version 
of  the  RTG  model,  the  method  of  solving  the  equations  of 
motion  of  a structure  with  viscoelastic  components  described 
by  RT  models  will  be  seen  to  be  a special  case  of  the 
following  solution  technique. 

The  reason  for  developing  a special  solution  technique 
for  the  equations  of  motion 

s 2[M]  {X  (s ) > + [K (s ) ] {X  (s ) } = (F  Cs ) } 6.7 

is  that  the  normal  method  of  decoupling  the  equations  of 
motion,  using  modal  analysis  to  construct  an  orthogonal 
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transformation  that  diagonalizes  the  mass  and  stiffness 
matrices,  is  not  applicable  because  the  stiffness  matrix 
of  the  structure,  [K(s)]  , contains  terms  that  are  depen- 

dent on  the  Laplace  parameter,  s 

The  method  of  solution  for  Eqn  6.7  is  an  extension  of 
the  method  proposed  by  Foss  to  decouple  the  equations  of 
motion  for  a structure  with  non-proportional  viscous 
damping  (Ref  26). 

[M]  { x (t ) } + [C]  {x(t)>  + [K]  { x (t ) } = { f (t) } 7.1 

Non-proportional  damping  occurs  when  the  damping  matrix  [C] 
is  not  a linear  combination  of  the  mass  and  stiffness  matrices 
of  the  structure. 

[C]  f aj[M]  + a2[K]  7.2 

At  present  there  is  no  general  method  of  constructing  an 
orthogonal  transformation  for  three,  real,  square,  symmetric 
matrices  when  each  of  the  matrices  is  not  a linear  combination 
of  the  remaining  two.  Consequently,  Foss  posed  the  equations 
of  motion  for  non-proportional  damping  in  terms  of  two  real, 


symmetric  matrices. 


The  lower  set  of  the  partitioned  matrix  equations  is  the 
equations  of  motion  of  the  structure  and  the  upper  set  of 
matrix  equations  is  satisfied  identically.  The  equations  of 
motion  as  posed  in  Eqn  7.3  are  readily  decoupled  and  solved. 

To  solve  the  equations  of  motion  of  the  structure  con- 
taining elastic  and  viscoelastic  components,  Eqn  6.7,  the 
equations  are  posed  in  terms  of  two  real,  square,  symmetric 
matrices.  To  begin,  one  multiplies  the  equations  of  motion 
by  each,  distinct  term  appearing  in  the  denominators  of  the 
elements  of  the  stiffness  matrix,  [K(s)]  , 

[D(s)s2[M]  + D(s)[K(s)]]  {X  (s ) } = D(s)  {F(s)}  7 


where 

K p 

N n B , N n 6 

DCs)  = n (1  + l a s nk)  . n (1  + l b s-"P)  7 

n=l  k=l  nk  n=l  p=l  np 

assuming  that  there  are  N different  viscoelastic  materials 
in  the  structure.  Multiplying  D(s)  times  [K(s)]  , the 

stiffness  matrix,  produces  a matrix,  [Kp(s)]  , that  has 

no  terms  in  s appearing  in  the  denominators  of  its  elements 

D(s)  [ K (s ) ] = [ Kp (s ) ] 7 

In  fact,  all  of  the  elements  of  [KD(s)]  are  constant  terms 
plus  terms  containing  s raised  to  real,  positive  powers. 
Also  note  that  the  matrix  D(s)s2[M]  has  elements  which 


are  suras  of  terms  containing  s raised  to  real,  positive 
powers . 


The  equations  of  motion  are  now  expressed  as 

[ Z (s ) ] {X  (s ) } = D (s) {F  (s)  } 7.7 

where 

[Z  (s ) ] = [D(s)s2[Mj  + CKD(s)]]  7.8 

At  this  point  in  the  development  the  real,  positive  exponents 
of  s appearing  in  the  matrix  [Z(s)]  are  taken  to  be 
rational  as  well.  Had  any  of  the  exponents  been  initially 
irrational,  they  are  replaced  by  their  rational  approximations 
to  as  many  significant  digits  as  desired.  Since  all  of  the 
exponents  in  [Z(s)]  are  rational,  the  matrix  may  be 
expressed  as 

[Z  (s)]  = [ [M]  l c.s3/m  + \ [KjsVm  ] 7.9 

j = 2nr  1=1 


where  m is  the  smallest  common  denominator  of  the  exponents 
of  s in  [Z(s)]  and 


[KD(s)]  = l [K  ]s£/m  7.11 

u £ = 0 1 

where  [K  ] is  symmetric  and  some  of  the  [K  ] and  c. 

X.  I J 

appearing  above  may  be  zero. 

Using  Eqn  7.9,  the  equations  of  motion  of  the  structure 
become 


[ [M]  l c.sj/m  + l [K,]s£/m]  ( X C s ) } = D (s ) {F  (s ) } 7.12 


j = 2m 


J 


i=0 


and  expressed  in  terms  of  one  index  of  summation,  they  are 


J j / 

I [ [M]c.  + [K  ]]s  m {X  (s ) } = D(s) {F (s) } 

j = 0 J J 


7.13 


or 


/m 


J 

l [A,]s  'm  {X  (s ) } 
j = 0 J 


= D (s  ) {F  (s ) } 


7.14 


where 


IXj]  = [ [M]c.  + [K.]  ] 


7.15 


and  again  recognizing  that  some  of  the  Cj  and  [K^]  are 
zero. 

The  equations  of  motion  as  given  in  Eqn  7.14  are  now 
posed  in  terms  of  two  real,  square,  symmetric  matrices. 
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7.16 

Note  that  each  matrix  [Aj]  , j = 0 ,1 , 2 , . . . , J , is  real, 
square  and  symmetric  because  they  are  linear  combinations 
of  [M]  and  [Kj]  . It  follows  that  the  two  matrices 
containing  [A^ ] in  Eqn  7.16  are  also  real,  square. and 
symmetric.  Also  notice  that  the  lowest  set  of  partitioned 
matrix  equations  in  Eqn  7.16  are  the  equations  of  motion  of 
the  structure  as  given  by  Eqn  7.14,  and  that  all  of  the  upper 
sets  of  matrix  equations  in  Eqn  7.16  are  satisfied  identically. 

The  equations  of  motion  as  posed  in  Eqn  7.16,  referred 
to  as  the  expanded  equations  of  motion,  can  be  decoupled 
using  an  orthogonal  transformation.  The  general  form  of  the 
expanded  equations  of  motion  is 


{0.} 

(0.} 

(0.) 


(0.) 

{0.} 

D(s){F(s)} 


f 


which  is  Eqn  7.16  expressed  in  more  compact  notation.  The 
orthogonal  transformation  is  constructed  from  the  eigen- 
vectors associated  with  the  eigenvalue  problem  for  the 
expanded  equations  of  motion. 

Xn[M]Un}  + [K]Un}  = {0.}  7.18 

The  eigenvectors  {<j>n>  are  used  to  construct  the  orthogonal 
transformation  matrix  [$]  in  the  normal  fashion,  and  the 
resulting  transformation  of  variables  is 

{X  (s ) } = [i]{a(s)}  7.19 

Substituting  this  transformation  into  Eqn  7.17  and  premulti- 
plying  the  equation  by  [<!>]  produces 

j 

s /m[i]T[M][i]{aCsn  + |>]T[K][>]  a(s)  = [i]T{F  Cs) } 7.20 

To  demonstrate  that  Eqn  7.20  is,  in  fact,  the  decoupled 
expanded  equations  of  motion,  one  uses  the  fact  that  eigen- 
vectors of  the  expanded  equations  of  motion  are  orthogonal 
with  respect  to  [M]  and  [K] 


Wj}T[M]Un} 


0.  j { n 


7.21 

7.22 


Equation  7.20  then  reduces  to 


1/ 

s m[^m-^]{a(s) } + Ok--J{a(s)}  = [*]1{F(s)} 


' -.T  , 


n 


n 


7.23 


where  ] and  [~~-k^]  are  diagonal  matrices  of  the 

modal  constants  m and  k , respectively. 

n n r 


% “ <*„»  M<  n> 


7.24 


k„  - Un)T[K]{in} 


7.25 


Premultiplying  Eqn  7.23  by  * or  equivalently 

yields 


m. 


n 


(a(s)  } + ["— -Of  a (s)  } 
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hm-J[j]T(F(s)) 
1 n 


7.26 


The  r".tio  of  the  nt^1  modal  parameters,  k /m  , is 

r ’ n n 

minus  the  n1  eigenvalue  of  the  expanded  equations  of  motion. 


Vmn  ' A 


7.27 


Premultiplying  Eqn  7.18  by  { ^ n } produces 


y*n}TM{*n}  + Un}T[^]Un}  = 0 


7.28 


X m + k =0 
n n n 


7.29 
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from  which  Eqn  7.27  follows.  As  a result,  the  equations  of 
motion  can  be  further  simplified  to 


1 / 

s { a (s ) } 


C — 3 { a (s  ) } 


["SpJ[$]T{F(s)}  7.30 

n 


From  Eqn  7.30  it  follows  that  the  expression  for  the 
Laplace  transform  of  the  nt^'  modal  coefficient,  an(s)  > is 


an(s)  * 


Un>T(FCs)> 
mnfs  /"'-xn) 


7.31 


This  expression  for  the  modal  coefficient  can  be  further 
simplified  by  noting  the  general  form  of  the  eigenvector 
{ <j>n } associated  with  the  expanded  equations  of  motion 


r 


;£'1(v 


x^“2{ 4>  } 
n Tn 


(*„) 


Xn<V 


Xn{V 


{V 


J 


7.32 


where  {<}>„}  and  Xn  are  the  solutions  of  the  eigenvalue 
problem  associated  with  the  original  equations  of  motion  in 
the  form  of  Eqn  7.14. 


[ l [A.]xj[]{*}  = {0.}  n = 1,2,3, ...,N  7.33 

j = 0 ->  n n 


The  general  form  of  the  nt^1  eigenvector,  { 4-  } , can  be 

verified  by  direct  substitution  in  Eqn  7.18  when  [M]  and 

[K]  are  expressed  in  terms  of  the  matrices  [Aj ] as 

indicated  in  Eqn  7.16.  The  lowest  set  of  resulting, 

partitioned  matrix  equations  produces  Eqn  7.33,  and  the 

upper  sets  of  the  partitioned  matrix  equations  a:e  satisfied 

identically.  Consequently,  the  numerator  of  the  n*^  modal 

T ' 

coefficient,  { <i>n } {F(s)}  , is 


Un}T{F(s)  } 


r ■> 

T 

r ^ 

1 

(0.  } 

■J-  2 , 

xn  {*n} 

{0.  } 

X^'3{4  } 

n n 

(0.  } 

'V*  , *\j  *\j  . *\, 

'Vf  • *\/  ^ . Or 


Xn{*n} 

W 

{*n> 


{0.  } 

{0..} 


J 


^ D(s){F(s)}  J 
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or 


(in}T{F(s)}  = Un)T{F(s)}D(s) 


and  the  n1"*1  modal  coefficient  reduces  to 


an(s) 


Un>  {F(s)}D(s) 


m (s  ) 

n v n' 


The  Laplace  transform  of  the  displacement  response  of  the 
structures  follows  from  Eqn  7.19  and  takes  the  form 


N 

(X  (s ) } = l an(s)U  } 

n=  1 n n 


or 


N U_)J{F(s)}D(s) 
(X  (s ) } = l — — j- 

n=1  m (s  /m-X  ) 
n v n 


{^n> 


where  N is  the  order  of  the  matrices  [M]  and  [K]  in 
the  expanded  equations  of  motion. 

The  order  of  the  expanded  equations,  N , can  he 
very  large.  From  Eqn  7.16  it  is  clear  that  the  order  of 
the  expanded  equations  of  motion  is  equal  to  5 , the 

order  of  the  matrices  [A^]  , times  J where  J is 

defined  in  Eqn  7.10. 


7.35 


7.36 


7.37 


7.38 
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N 


6 • J 


7.39 


From  Eqn  7.10,  it  is  clear  that  is  the  largest 

exponent  in  the  expression  s2D(s)  which  is  2 + g , 
where  g is  the  largest  exponent  in  D(s)  . Therefore, 

J = m(2+g)  7.40 

and  the  order  of  the  expanded  equations  is  seen  to  be 

N = 6m (2  + g)  7.41 

Note  that  if  m , the  smallest  common  denominator  of  the 
rational  exponents  of  s , in  the  original  equations  of 
motion  is  large,  the  order  of  the  expanded  equations  is 
quite  large  for  a structure  with  anything  more  than  a very 
modest  number  of  degrees  of  freedom.  However,  the  solutions 
to  the  expanded  equations  of  motion  can  be  obtained  by  using 
numerical  methods  that  do  not  involve  the  manipulation  of 
the  expanded  equations,  as  will  be  shown  in  the  following 
section . 

Before  proceeding,  it  should  be  pointed  out  that  the 
equations  of  motion  of  a structure  containing  both  elastic 
and  viscoelastic  components  can  be  solved  given  that  a 
finite  element  formulation  of  the  equations  of  motion  is 
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possible  and  that  an  RT  or  RTG  model  exists  for  each 
viscoelastic  material  in  the  structure.  The  general  form 
of  the  solution  technique  for  the  equations  of  motion 
containing  only  RT  models  for  the  viscoelastic  components 
is  identical  to  the  above  development  except  that  D(s) 
is  set  equal  to  one. 


VIII.  Calculating  the  Laplace  Transform 
of  the  Structural  Response 


At  this  point,  it  is  clear  that  solutions  to  the  equations 
of  motion  for  the  structure  containing  elastic  and  visco- 
elastic components 


N {*  } (F(s)}D(s) 
(X(s)}  = l -.5— 1 

n=1  m m~X 

mn  (s  n 


<*n> 


7.38 


are  difficult  to  calculate  from  the  expanded  equations  of 
motion,  because  of  the  large  order  of  the  matrices  in  the 
expanded  equations.  As  a result,  an  alternative  method  of 
obtaining  the  solutions  is  required  if  the  finite  element 
formulation  of  the  equations  of  motion  using  generalized 
derivative  models  is  to  be  a useful  tool  to  the  engineer. 

The  alternative  method  adopted  here  is  a combination 
of  iterative  schemes  used  to  obtain  the  eigenvalues, 

Xn  , and  eigenvectors,  { 4>n } , associated  with  the  original 

equations  of  motion. 

[ I [A.]xJ]{*J  = {0.}  , n = 1,2,3 N 7.33 

j=0  3 n n 


Recall  that  the  number  of  distinct  homogeneous  solutions  to 
the  equations  of  motion,  N , is  dependent  on  the  smallest 
common  denominator  of  the  exponents  in  the  equations  of 
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motion,  m ; the  largest  exponent  in  the  product  of  the 
denominators  terms  of  the  Laplace  transform  of  the  RTG 
models  in  the  equations  of  motion,  6 ; and  the  number 

of  degrees  of  freedom  of  the  structure,  6 

N = 6m(2+0)  = 2m6  + 8m6  7.41 

Using  an  iterative  scheme  based  on  the  homogeneous 
form  of  the  equations  of  motion  in  Eqn  6.7,  26m  of  the 

homogeneous  solutions  are  obtained. 

[xJm[M]  + [K(Xn)]Un}]  = {0.}  8.1 

The  iteration  for  the  solutions  centers  on  calculating 

successive  estimates  of  and  using  the  scheme 

n n & 

[-2m(P+n[M]  + [K(XiJP))]]{4>n}(P+n  = {0.}  8.2 

X^  is  the  p*"^1  estimate  of  X . ^2m(P+l) 

n r n n 

(p  + I)**1  estimate  of  X^m  and  {<)>n}^P+^  is  the 
(p  + l)**1  estimate  of  { <(in } 

Given  a value  of  X^P^  > the  (p  + l)**1  estimates 
of  x‘m  and  {4>n)  may  he  calculated  using  matrix  iteration 
or  any  other  method  that  is  appropriate  to  obtain  the 
solution  to  the  eigenvalue  problem 
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[ B ] { <X> } 


W { 4> } 


8.3 


where  B is  complex.  Note  that  Eqn  8.2  can  be  expressed  as 


[K(X<P))]'1[M]Un}(P+1) 


1 i ± i(P+l) 

{V 


n 


8.4 


which  is  of  the  same  general  form  as  Eqn  8.3. 

The  iterative  scheme  as  it  appears  in  Eqn  8.4  is  only 

useful  in  obtaining  the  eigenvalue  with  the  smallest 

magnitude,  \ ^ , and  the  associated  eigenvector,  { 4>  ^ } 

To  obtain  the  other  eigenvalues  and  eigenvectors,  the 

iteration  scheme  is  modified  to  allow  the  scheme  to  converge 

on  the  larger  eigenvalues  and  associated  eigenvectors. 

ti  h 

For  instance,  the  scheme  used  to  obtain  the  L eigen- 
value, , and  the  Lth  eigenvector,  {$L)  , is  (Ref  27) 


[Kti'1’’)]-1  * I 
L 1=1 


1 {^l}(P){^SL}(Py 


2m  (Pi 


[M]Wl> 


(P+1) 


-2m (P+1) 
AL 


<*L} 


(P+1) 


8.5 


The  terms  in  the  summation  on  the  index  I subtract  the 
components  of  the  first  L - 1 eigenvectors,  associated 
with  eigenvalue  problem, 


'A  ■ ' 
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r 


1 


[K(x[P))]_1[M]{^}(P)  = - U£)(P)*-1,2,3 L-l 

A } 

* 8.6 

from  the  successive  approximations  of  { <j> L } ^P+ 1 ^ calculated 
when  matrix  iteration  is  used  to  solve  Eqn  8.5. 

Given  that  matrix  iteration  has  successfully  produced 
x[m  and  (<J)jj}^P+^^  , one  can  take  advantage  of  the  orthogon- 

ality of  the  mode  shapes 


[M]{*l}(P+1)  = 0.  1-1,2,. ..,L-1 


8.7 


to  demonstrate  that 


CP) 


j1^{V(f)Iw('L»(,*1) 


and  Eqn  8.5  reduces  to 


{0.  } 


8.8 


[KCx'P))]'1[M]{*L)(Ptl) 


1 ,CP*1) 

:2m (P+1)  {*L} 

AL 


8.9 


I 


4 


> 

I 


! 


r 


Thus , {<J>L>^P  + 1^  and  A2m(P+l)  are  fact  the  (p  + l)t^ 

approximation  of  the  solution  to  the  equations  of  motion. 

On  the  other  hand,  { ^ ^ P and  a£P^  » 1 = 1 * 2 L-l 

are  not  in  any  sense  approximations  of  the  solutions  to 

the  equations  of  motion.  However,  they  are  the  first  L-l 

eigenvalues  and  eigenvectors  associated  with  Eqn  8.9.  This 

fP) 

can  be  seen  by  comparing  Eqn  8.6  with  Eqn  8.9.  and 

f PI 

A;  1 can  be  calculated  using  matrix  iteration. 
l ° 


££?**■  - iw* 


i 

I 


Notice  that  in  continuing  the  iterative  processes  in 

Eqn  8.2  or  8.5  that  the  (P  + 2)*^  approximations  of  the 

solution,  ^2m(P+2)  an(j  ^ j(P+2)  are  based  on  the 

n 71  ’ 


rn 


1/ 


- rp+i  ■»  J 

value  of  J However,  the  function  Z has  2m 


2m 


branches.  Given  a value  of  ^2m(P+l)  one  can  caicuia^e 

n ’ 

2m  values  of  X^P+1^  , one  value  for  each  branch  of 


2m 


So,  when  using  Eqns  8.2  and  8.5  to  obtain  solutions 
of  the  equations  of  motion,  it  is  necessary  to  choose  one 


1/ 


2m 


branch  of  Z “m  to  calculate  the  6 eigenvalues  and  6 
eigenvectors.  Then  another  branch  of  Z m is  chosen  and 
6 other  solutions  are  obtained.  This  process  is  continued 
until  all  2m  branches  have  each  been  used  to  calculate 
6 solutions  producing  a total  of  2m5  homogeneous  solutions 
to  the  equations  of  motion.  The  general  form  of  the  equations 
that  has  these  2m6  homogeneous  solutions  is 


[X^k)[M]  + [K(X.fvO]3Uwv'|}  = {0*}  j 1,2 6 


jOO 


pj  00 


k=l , 2 , . . . , 2m 


where  the  subscript  k denotes  the  branch  of  Z 
which  the  relation 


1/ 


2m  on 


8.10 


f \ 2m  11/2m  = ~ 

1 j(k)J  Xj(k) 


8.11 


is  valid. 
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The  remaining  gm,$  of  the  N homogeneous  solutions 
to  the  equations  of  motion  are  determined  using  Eqn  7.41, 
an  iterative  scheme  based  on  the  homogeneous  form  of  Eqn  7.12.* 


[[M]  1 c.V  * I [K  ]xj]{6  } 

j = 2m  J n o = n *■  « n 


{0.} 


8.12 


1= C 


The  iterative  scheme  is 


[M  (CjXn 


CP+lJrJ-l 
An 


(P)  J-l.  ..(P) 

+ y xc.xJ  ) 

. L~  in 
j = 2m  J 


£ = 0 


n 


{0.  }8. 13 


or,  in  more  compact  form 


[ [M]  Q (X^P+1  ^ x^P^ ) + [KD(X15P))]]{^}(P+1)  = {0.}  8.14 


where,  as  before,  X J-P^  is  the  P*^  estimate  of  the  eigen- 
value, ^P  + estimate  of  the  eigen- 

value and  {41  } ^P+^  is  the  (P+l)1"*1  estimate  of  the 
n 

- r p ") 

eigenvector.  Given  a value  for  x^  ' , one  can  calculate 

rp+n  - rp+1 1 * rpi  . . . 

{d>  } 1 J and  Q(XV  1 , \K  ')  using  matrix  iteration 

n ^ n ’ n 0 

or  any  other  method  suitable  for  the  solution  of  Eqn  8.3. 


* Note  that,  if  the  equations  of  motion  contain  only  RT  models  for 
the  viscoelastic  materials,  that  D(s)  is  one  and  e,  the  largest 
exponent  of  s in  D(s),  is  zero.  Hence,  the  total  number  of  homogen- 
eous solution,  N,  is  2m6  as  seen  by  Eqn  7.41.  Since  Eqns  8.2  and 
8.5  provided  2m6  solutions,  one  can  return  to  Eqn  7.38  and  calculate 
structural  responses. 


I 
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Equation  8.14,  expressed  in  the  form  of  Eqn  8.3, 


is 


[KD(*J[p))r1[M]wn}(p+1) 


QCx 


1 r i ,(P+1)  o i r 

: (p+1  )'  ;'CH,  {V  8*15 


n 


’ xn  J 


f p + 1 'j  CP') 

At  this  point,  the  numerical  value  of  Q(XV  ,x 
r ’ ^ n n 

from  Eqn  8.25  is  used  to  calculate  x^P+^  with  the  relation 


Od(Ptl)Am) 

n * n ’ 


J n n 


(P) 


l c.J<p> 
j = 2m  ^ n 


8.16 


which  follows  from  Eqns  8.13  and  8.14.  However,  this  method 

of  calculating  x^P+^  assumes  that  are  on 

the  principal  branch  of  (x^m)  ^ • The  method  of 

calculating  X ^P  + ^ assuming  that  x ^P+^  and  X ^P^  are 

b n b n n 

both  on  the  k*^1  branch  of  (x6m)  is 

n 


’orx(p+1)  ).  !c  x(p) 

^ n (k)  ’Xn(k)J  i^1CjXn(k) 

CPTTJ 


c l3'1 
CjAn(k) 


1/ 


6m 


6m 


- ;(p+1) 
" n(k) 


8.17 


where  the  kt^1  branch  of  Z is  used  to  calculate  X^,+  ^ . 

n (k ) 

The  resulting  form  of  the  iteration  process  is 


CKD(Xn(k))]  1[M]{<f>n(k)}(P  ^ 


^f'(P  + l).  (P)'  Un(k)} 
QUn(k),Xn(k)J 


(P+1) 


8.18 
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* % 

- ~ •(>#  - 


where  successive  estimates  of  X 


n(k) 


are  calculated  using 


(k)  ’ ^ from  ^cln  8.18  and  then  using  Eqn  8.17  to 


calculate  x 


(P+1) 

n(k) 


When  using  matrix  iteration  to  solve 


for  Q(xn(k)  ^nCk)-1  and  { ^n  (k) } ^ in  Eqn  8-18>  one 

usually  obtains  only  the  Q with  the  smallest  magnitude 

and  its  associated  eigenvector,  • 

To  obtain  successive  estimates  of  the  other  Q’s 

with  larger  magnitudes  and  their  associated  eigenvectors, 

the  iteration  scheme  is  modified  as  before.  The  iterative 

• fp+ivrpi 

scheme  that  yields  Q (A  l (k)  ’ *L  (k)  ^ and  ^L(k)*  1S 


^KD(Xn(k)^ 


K1  {nm}(p){n(k)} 

£ii  o r a P+1 ^ A(p}  l 

a 1 qlA£(k)  ,At(k)J 


[ M ] { 4> . 


(P+1) 


( t ri,  A ) 


“"(p+i)'(p)  : l,pL(k) 

^UL(k)  ’XL(k)J 


(P+1) 


8.19 


(P+1 1 (PI  (PI 

where  Q(A  ^ (k)^  and  ^£(k)^  are  t'^le  s°lut;ions  t0 


£KD(xL(k))]  ^^Uk)  1 ^ 


o ("Av+bj*r){*i(u}(P)8-20 

QlAi(k)  ,A£(k)  J 


having  the  L - 1 smallest  values  of  Q . The  values  of 
^A£Pk)  ’ At(k)  ^ and  ^£(k)^P^  can  Ee  oEtafnec^  using  matrix 
iteration. 

Using  Eqns  8.17,  8.18,  and  8.19,  66m  homogeneous  solution 
to  the  equations  of  motion  are  obtained  for  each  branch  of 


* > - s 

. . iFiy 


1/  1/ 

Z . Since  there  are  pm  branches  of  Z Bm  , the 

iteration  process  should  produce  Bmfi  homogeneous  solutions. 

These  homogeneous  solutions  satisfy  the  homogeneous  equations 

of  motion. 


J 

[ [M]  X cix 

j = 2m  J 


POO 


[Kjx 


] { 4>  v 


tV  * P<k>  P<k> 


(0.  } 


P = 1 > 2 , . . . , 6 
k=l , 2 , . . . , 6m 
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The  two  iterative  schemes  used  to  obtain  the  N solutions 

to  the  homogeneous  equations  are  considered  to  have  converged 

~ f P+1 1 

when  successive  approximations  of  xn(k)^n(k)  an<*  *n(k)  ^ 
are  approximately  the  same  complex  number. 


: (p+1) 

Xn(k) 


1.0 


8.22 


The  general  convergence  criteria  of  the  iterative  schemes  are 
considered  beyond  the  scope  of  this  investigation;  however, 
the  schemes  have  converged  for  numerous  structures  considered 
by  the  author.  Homogeneous  solutions  for  an  example  problem, 
calculated  using  the  iterative  scheme  given  above  appear  in 
Table  8.1. 

Note  that  all  of  the  parameters  of  the  Laplace  transform 
of  the  structure’s  response  are  determined,  except  the  modal 
constant  m 


defined  by  F.qn  7.24. 


Table  8.1 

HOMOGENEOUS  SOLUTIONS  OBTAINED  USING 


THE  PROPOSED  ITERATIVE  SCHEMES  FOR  AN  EXAMPLE  PROBLEM 


.3300  .0825 

.0825  .3300 


[K(s)] 


4.0  + 


0.2  + 0,2s 

1.0  + 0.1s' 


-2.0 


- 2.0 


4.0  + 


0.01  + 0.01s 
1.0  + 0.1s*5 


Using  Eqns  8.2  and  8.5 


X1(1^  = 1.071768  + i 1.099794 


f 1.000  + i 0.0  1 

l 1.040  + i O.OI65J 


ll(2) 


ll(4) 


1(5) 


2 (1 ) 


2(2) 


1.073498  + i 1.028611  U 


fl.OOO  + i 0.0 
} =1 

Ll.001  + i 0. 


1(2) 


1.073498  - i 1.028611  U1(-3)) 


1.071768  - i 1.099794  ^ (4 ) > 


1.581998  + i 1.601989  ($2(1)} 


■1.584631  + i 1.547150  U2(2)} 


0238' 


fl.OOO  + i 0.0  | 
ll.OOl  - i 0.0238-1 
fl.OOO  + i 0.0  1 
1 1 . 040  - i 0.0165 / 


f 1. 000  + i 0.0  ) 
1-.9713+  i 0.0120-* 
fl.OOO  + i 0.0  ) 
1-1.002  + i 0.0240) 
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TABLE  8.1 
(Continued) 


9.99385  + i 0.0 


X2(4)  = -1.584631  - i 1.547150  U2(3)> 


0.03  - i 0.0240 


2 (5) 


1.581998  - i 1.601989 


U2(4)} 


| 1.000  + i 0.0 
l-  .971  - i 0.0120 


Using  Eqns  8.18,  8.19,  and  8.20 

Xi(3)  = -9.997418  + i 0.0  {<frl(3)} 


*2(3) 


{4>2(3)} 


7.24 


mn  = {*n}  [M]Un} 


The  eigenvector  of  the  expanded  equations  of  motion,  {<f>n}  » 

can  be  constructed  from  the  n**1  eigenvalue  and  associated 
eigenvector  of  the  original  equations  of  motion,  *n  and 
{ <(>n } , using  the  relation  given  in  Eqn  7.32.  This, 

coupled  with  the  general  form  of  [M]  given  in  Eqn  7.16, 
produces  an  expression  for  mn  which  takes  the  form 


mn  = {V  Cj=V*An 


j-1 


[Aj  33 ( ) 


8.23 


The  modal  constant  of  the  expanded  equations  of  motion, 
mn  , can  be  calculated  without  manipulating  the  expanded 
equations  of  motion. 

In  conclusion,  all  of  the  parameters  in  the  general 
form  of  the  Laplace  transform  of  the  structure's  response 
can  be  calculated  without  manipulating  the  expanded  equations 
of  motion,  given  that  the  iterative  schemes  outlined  above 
converge . 
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IX.  The  Existence  of  the  Structural 


Response  to  Impulsive  Loading 

Of  particular  interest  at  this  juncture  is  whether  or 
not  the  inverse  transform  of  the  Laplace  transform  of  the 
structure's  displacement  response  for  impulsive  loading 
exists.  In  fact,  the  inverse  transform  always  exists  and 
is  real,  continuous  and  causal. 

To  demonstrate  this,  one  starts  with  the  general  form 
of  the  transform  of  the  structural  response. 

7.38 

at  the  structure's 
applied  forces  is 

9.1 

The  Laplace 
is 

(F  (s ) } = {1.}  9.2 


{x  (s) ) = I <♦„  |T|r(s)tP(s) 


n=l 


mn(s"r7”'Xn) 


For  simultaneous  unit,  impulsive  loading 
degrees  of  freedom,  the  column  vector  of 

{f(t)}  = 6 (t  ) { 1 . } 


where  {1.}  is  a column  vector  of  ones, 
transform  of  the  column  vector  of  forces 
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-\\r-  ■; 

- 


and  the  transform  of  the  response  for  the  impulsive  loading 


tics))  - I <vTu->nts) 


nil  1/ ■ - «♦„> 

m (s  ra-X  ) 


The  inverse  transform  of  this  expression  always  exists, 
which  follows  from  a theorem  on  the  existence  of  the  inverse 
transform  (Ref  28).  Paraphrasing  the  theorem  in  terms  of 
the  notation  used  above,  it  states  that  the  inverse  transform 
of  (X(s)}  exists  and  is  real,  continuous  and  causal  when 


1.  (X(s)}  is  analytic  for  Re[s]>0  , 


2.  (X(s)}  is  real  for  s real  and  positive,  and 


3.  (X(s)}  is  order  s’Y  , where  y>1  , for  s 

large  in  the  right  half  s plane. 


(X(s)}  is  analytic  for  Re[s)>0  when  the  branch  cut 

1/ 

of  s 171  is  chosen  to  lie  along  the  negative,  real  axis 

A 

in  the  s plane  and  the  poles  of  (X(s)}  , which  occur  at 


and  do  not  appear  in  the  right  half  s plane. 


* A pole  in  the  right  half  s plane  indicates  that  a RT  or  RTG  model 
in  the  equations  of  motion  characterizes  the  viscoelastic  material  as 
generating  energy  instead  of  dissipating  energy. 


A * • • v 
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{X(s)J  is  real  for  s real  and  positive.  The  only 
quantities  appearing  in  Eqn  9.3  that  can  be  complex  are 

l/m 

{ } , A and  m , because  DCs)  and  s are  real 

n n n 

for  s real  and  positive.  When  { 4>n } , and  mR 

are  complex,  they  occur  in  conjugate  pairs.  Note  if  An 
and  are  a homogeneous  solution  to  the  equations 

of  motion 


J 

[ l 

j = l 


[Aj]^<*n} 


{0.  } 


7.33 


where  and  C <J>n ) are  complex,  that  their  conjugates  are 

homogeneous  solutions  to  the  equations  of  motion. 


J 

C l 
j = l 


{0.  } 


9.5 


Equation  9.5  follows  directly  from  the  complex  conjugate  of 
Eqn  7.33.  Since  homogeneous  solutions  occur  in  conjugate 
pairs,  the  modal  constant  mn  occurs  in  conjugate  pairs. 


m. 


n 


c $ } 1 [ y j • x 
n Lji1J  n 


3-1 


[Aj]]Un} 


8.23 


Un}  Cj^1j,Xn 


j-1 


CAj  ]3  {4»n> 


9.6 


It  follows  directly  that  when  the  terms  in  Eqn  9.3  are 
complex,  they  occur  in  complex  conjugate  pairs  and  (X(s)} 
is  real  for  s real  and  positive. 


oSSfe* 


69 


{X(s)}  also  satisfies  the  third  and  last  condition 

A 

placed  on  the  transform.  To  show  that  (X(s)}  is  order 
_2 

s for  s large  in  the  right  half  s plane,  one  uses 
the  transformed  equations  of  motion  as  they  appear  in 
Eqn  6.7. 


s*[M]{X(s)}  + [K(s)]{X(s)}  = { F (s ) } 


6.7 


The  transformed  equations  of  motion  for  simultaneous,  unit 
impulsive  laoding  is 


[s*[M]  + [K(s)]]{X(s) } = {1.} 


9.7 


The  only  terms  in  [K(s)]  , other  than  the  constant  terms, 

are  those  terms  proportional  to  p*(s)  and  x*(s)'  from 
the  RT  and/or  RTG  models  of  the  viscoelastic  materials.  The 
general  forms  of  p*(s)  and  X*(s)  appear  in  Eqn  4.29  for 
the  RTG  model  and  Eqn  4.30  for  the  RT  model.  As  a direct 
result  of  the  general  forms  of  p*(s)  and  X*(s)  , 

Eqn  9.7  reduces  to 


s^[M]{X(s)} 


(1.  } 


9.8 


for  s large.  Since  the  elements  in  the  mass  matrix  are 

2 


constant,  {X(s)}  must  be  order  s 

-2 


(X(s)}  is  order  s 
plane . 


Therefore , 
for  s large  in  the  right  half  s 
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Having  now  established  that  response  of  the  structure 
to  impulsive  loading,  (x(t)}  , 

{x(t)}  = L* 1 {X  (s ) } 9.9 

exists  and  is  real,  continuous  and  causal,  the  next  issue 
to  be  addressed  is  the  calculation  of  the  inverse  transform, 
the  topic  of  the  next  section. 


I 


X.  The  Calculation  of  the  Response  to  Impulsive  Loading 


The  final  step  in  determining  the  impulse  response  of 
the  structure  is  to  calculate  the  inverse  transform  of 
the  Laplace  transform  of  the  response  to  impulsive  loading. 


(X  (t ) } = L- 1 {X (s ) } 


■1 


N 

l 

n=l 


}T{1.  >D(s) 

• P-17  — r Un} 

mn(s  m-Xn 


10.1 


The  inverse  transform  integral 


L-1{X(s)} 


1 ,Y+1°°  st  » 

— J est{X  (s ) } ds 

2ni  y - i“ 


10.2 


is  evaluated  using  the  residue  theorem  from  the  calculus 
of  a complex  variable. 

The  closed  contour  of  integration,  used  in  conjunction 

with  the  residue  theorem,  is  given  in  Fig  10.1.  The  contour 

is  divided  into  six  segments  and  the  direction  of  integration 

along  each  segment  is  indicated  by  the  arrows.  Segments  3, 

4,  and  5 of  the  contour  are  required,  because  the  branch  cut 

1/ 

and  branch  point  of  the  function  s m are  taken  to  be 
along  the  negative  real  axis  and  at  the  origin  of  the  s 
plane,  respectively. 


- 
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The  residue  theorem  states  that  the  integral  along  the 


closed  contour,  divided  by  2ui  , is  equal  to  the  sum 
of  the  residues  of  the  poles  of  the  integrand.  In  this 
case,  the  statement  of  the  residue  theorem  translates  into 


/ {X(s)}eSt  ds  = - yL-  l ! {X  (s)  }est  ds  + jV  10.3 

1 k=2  k j 

where  the  circled  index  indicates  the  contour  over  which  the 
integration  is  to  be  performed,  and  b^  are  the  residues 
of  the  poles  of  {X(s)}  enclosed  by  the  contour. 

The  integrals  on  the  segments  of  the  closed  contour 
are  evaluated  for  the  case  where  the  length  of  segment  1 
is  extended  indefinitely  in  the  positive  and  negative 
imaginary  directions. 


I 


showing  that  one  need  evaluate  the  right  sid?  of  Eqn  10.5 
to  obtain  the  inverse  transform. 


To  maintain  the  continuity  of  the  closed  contour,  the 
radii  of  segments  2 and  6 are  increased  indefinitely,  and 


segments  3 and  5 are  extended  indefinitely  in  the  negative 
real  direction.  When  the  radii  of  contours  2 and  6 are 
increased  indefinitely,  it  can  be  demonstrated  that  the 
resulting  value  of  the  integrals  along  these  two  segments 
is  zero.  This  follows  directly  from  the  fact  that,  for 

a ^ o a.  y. 

large  s , (X(s)}  is  order  s . Hence,  the  i 

~ _ 2 
component  of  (X(s)}  is  order  s 


Xt(s) 


|s|  >>  1 


Therefore,  the  asymptotic  expression  the  integral  on 
segment  2 as  its  radius  increases  is 


/ X£( s)eStds 
2 


Re 

/ 1, 

Re 


1 IT 


ds 


10.6 


10.7 


Expressing  s in  polar  notation,  the  relation  is 

* (-<-  " Re^t 

/ Xt(s)e  ds  i J c8  Rie  de 

2 eo  R e 

I s I >>  1 


10.8 


The  magnitude  of  an  integral  is  less  than  or  equal  to  the 
maximum  value  of  the  magnitude  of  the  integrands,  M , 
times  the  length  of  the  path  of  integration,  L 
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1 1 


I 


4 


i 


,*  aRel9t  R ie.. 

O 


< M • L 


10.9 


The  maximum  value  of  the  magnitude  of  the  integral  is 

.Yt 


M = 


ive 


10.10 


and  the  length  of  the  path  of  integration  in  radians  is 


L = it  - 0 


10.11 


The  resulting  bound  for  the  magnitude  of  the  integral  is 

.Yt 


M • L = 


R 


which  is  clearly  zero  for  finite  time  in  the  limit  as  R 
becomes  large  or  equivalently,  in  the  limit  as  s becomes 
large.  Since  the  magnitude  of  the  integral  is  zero,  both 
the  real  and  imaginary  parts  of  the  integral  are  zero. 


lim  J X(s)estds 
R - 


= 0 . + iO  . 


The  expression  for  the  magnitude  of  the  integral  along 
segment  6 is 


i 
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-0o  eReiet  u 

I c - Re  ide  < M • L 

-u  1 „ 2 _ i 2 e 


10.14 


and  the  proof  that 


lim  / X (s)estds  = 0.  + iO, 


10.15 


follows  the  steps  given  in  Eqns  10.10  through  10.12. 

The  integral  along  segment  4 of  the  contour  is 
evaluated  for  the  case  where  the  radius  of  the  contour 


shrinks  to  zero.  The  asymptotic  expression  for  the  inte- 
grand, the  component  of  Eqn  9.3,  for  small  s is 


N (♦_}4{1.  } 

x , Cs ) ^ l —2 *_ 

n-l  mn(-xn) 


I S I <<  1 


10.16 


111  <<  I A. 


where 


|>nJl  is  the  ft1  component  of  { 4>n } , and 


DCs) 


10.17 


I s I <<  1 


Tf— 

s m-x 


n 


I s I <<  Unl 


10.18 


The  asymptotic  expression  for  the  integral  on  segment  4 is 


J X (s)estds  * / l — 

4 


N {6  }X{1.} 


4 n=1  mn(-Xn) 


*nteStds 


10.19 


I s | <<  1 
Is!  <<  |x. 


n 


In  polar  notation,  where 


s = pe 


ie 


10.20 


the  expression  the  integral  on  segment  4 is 


f X „(s)estds  * / , - s 

J J i n=l  m ( - X ) 

4 * 


N {4.  } { 1 . > ie.  . Q 
l -H-.- epe  tipel6de 


n' 


10.21 


Again,  the  magnitude  of  the  integral  on  segment  4 is  less  than 
or  equal  to  the  maximum  magnitude  of  the  integrand,  M , 
times  the  length  of  the  path  of  integration,  L 


<♦«>  Peietl 


I 1 


" n«l  mn(-xn) 


pe  t • ie, 
e ipe  de 


< M • L 


10.22 
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The  maximum  magnitude  of  the  integrand  is 


M 


N {*  }T{1.} 

I — 

n=1  v-v 


and  the  length  of  the  path  of  integration  is 


L = 2n 


10.23 


10.24 


The  resulting  expression  for  the  upper  bound  on  the  magnitude 
of  the  integral  on  segment  4 is 


/ X£(s)estds 
4 


N {*}T{1.} 

I — 

n=1  “nC-xn) 


ept2irp  10.25 


I s | <<  1 

I s I <<  Unl 


and  clearly  as  p goes  to  zero,  or  equivalently  as  s goes 
to  zero,  the  upper  bound  on  the  magnitude  of  the  integral 
goes  to  zero.  Therefore, 

lim  / X (s)estds  = 0.  + iO.  10.26 

p-o  4 

It  has  been  demonstrated  that  the  integrals  on 
segments  2,  4,  and  6 are  zero,  which  reduces  the  expression 
for  the  inverse  transform  given  in  Eqn  10.5  to 
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Y+i® 

— / (X  (s ) }estds 

2iri  y-i® 


— If  {X(s)}estds  + J>.  10.27 
!iri  k=3,5  r j ^ 


The  inverse  transform  is  seen  to  be  the  integrals  along  the 
branch  cuts,  segments  3 and  5,  plus  the  sum  of  the  residues 
of  the  poles,  bj 

The  expression  for  the  integral  along  segment  3 is 


/ X (s)estds  = 


D • i tt 

f v . i tt re  t itt, 

/ X (re  )e  e dr 

R 1 


10.28 


where 


re 


1 IT 


10.29 


This  expression  simplifies  to 


/ X (s)estds  = 


/ X (re11T)e'rtdr 
P * 


10. 30 


In  the  limit  as  the  radius  of  contour  2,  R , goes  to 
infinity  and  as  the  radius  of  contour  4,  p , goes  to  zero, 
the  integral  along  segment  3 becomes 


/ X (s)estds  = / X fre11T)e'rtdr 


lim 

p-»0 

R->-<» 


10.  31 


Similarly,  the  expression  for  the  integral  along  segment  5 
in  the  limit  as  the  radius  of  contour  6 goes  to  infinity 
and  the  radius  of  contour  4 goes  to  zero  is 


pio  J X fs)estds  = - / X (re‘in)e'rtdr  10.32 

R-"“  5 0 

The  sum  of  the  integrals  along  segment  3 and 
segment  5 is 


£ J X (s)estds  - 
k=3,5  k 1 


j"  (X  (reiir)  - X.(re'ilT))e*rtdr  10.33 
0 1 


Noting  that  X frelir)  and  X (re'1")  are  complex  conjugates 

jL  ** 

of  each  other,  the  expression  simplifies  to 


oo  . 

I J XA(s)estds  = - 2i  Im/  Xa(re‘lff)e‘rtdr 

k=3,5  k ° 

The  only  terms  in  the  inverse  transform,  Eqn  10.26, 
remaining  to  be  evaluated  are  the  residues  of  the  poles  of 
Xt(s)eSt  . By  definition,  the  residues  of  the  poles  of 
X4(s)est  are 


10.34 


or 


Before  taking  the  limit  to  evaluate  the  residue,  it  should 
be  pointed  out  that 


which  simplifies  the  expression  for  the  residue  to 


Evaluating  the  limit  produces 


Having  evaluated  the  residues  of  the  poles  of  the 
integrand  in  the  inversion  integral,  the  evaluation  of  the 
inverse  transform  is  complete.  The  general  expression  for 
the  displacement  response  of  the  structure  for  simultaneous, 
unit,  impulsive  loading  at  all  its  degrees  of  freedom  is 

x.(t)  = — Im  [/  X (re’11T)e’rtdr] 

IT  0 £ 

mX1!1"1  { <j>  • }T{1 . }D  (X^)  £t 

+ l —2 2 2—  * e 2 1.0 . 41 

j m.  ■> 

This  expression  is  based  on  Eqn  10.26,  using  Eqns  10.12, 

10.34,  and  10.40  for  the  appropriate  substitutions. 

Note  that  the  response  of  the  structure  has  two  parts, 
one  part  being  a sum  of  decaying  sinusoids  and  the  other 
part  an  integral  that  decreases  with  increasing  time.  The 
integral  does  not  decrease  exponentially,  because  it  is 
asymptotic  to  t’a  for  large  t , where  a is  greater 
than  one.  Therefore,  the  integral  dominates  the  response 
for  time  long  after  the  loading.  This  component  of  the 
response  describes  the  non-oscillatory  return  of  the  struc- 
{ ture  to  its  unloaded  equilibrium  position. 

In  summary,  the  response  of  the  structure  to  impulsive 
loading  can  be  calculated  using  contour  integration  to 
evaluate  the  inverse  transform.  Having  obtained  the  response 
to  impulsive  loading  as  a function  of  time  makes  possible 


* 


determination  of  the  response  of  the  structure  to  general 
loading  conditions  using  convolution.  In  essence,  the 
calculation  of  the  response  of  the  structure  is  no  longer 
tied  to  the  use  of  Laplace  or  Fourier  transforms.  In  other 
words,  the  response  of  the  structure  to  loading  time  histories 
for  which  transforms  do  not  exist  can  be  calculated  using 
the  response  to  impulsive  loading  and  the  convolution  integral. 


-1 
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XI . Summary  and  Conclusions 


The  generalized  derivative  of  fractional  order  is  a 
mathematical  operator  well-suited  for  describing  the 
frequency-dependent  mechanical  properties  of  viscoelastic 
materials.  As  shown  in  Section  V and  Appendix  C,  the 
generalized  derivative  constitutive  relations  for  the 
materials  are  capable  of  describing  the  frequency -dependent 
stiffness  and  damping  of  the  materials  over  several  decades 
of  frequency,  as  observed  under  conditions  of  sinusoidal 
motion. 

Moreover,  as  demonstrated  in  Section  V,  the  constitutive 
relation  for  3M-467  predicts  accurately  the  non-periodic 
(transient)  response  of  the  material  as  observed  in.  the 
laboratory.  In  addition,  the  generalized  derivative  model 
for  3M-467  performed  remarkably  better  than  the  Voigt  model. 

The  generalized  derivative  RT  and  RTG  viscoelastic 
models  enabled  the  formulation  of  the  viscostiffness  matrix 
of  the  viscoelastic  finite  element.  This  led  to  the 
successful  formation  of  the  equations  of  motion,  for  a 
structure  containing  both  elastic  and  viscoelastic  components, 
which  can  be  decoupled  and  solutions  obtained  using  modal 
analysis.  As  demonstrated  in  Section  VIII,  the  parameters 
of  the  solutions  can  be  determined  using  the  iterative 
numerical  schemes  presented. 


Finally,  it  was  demonstrated  that  the  response  of  the 
structure  to  impulsive  loading  always  exists  and  is  continuous, 
real  and  causal.  The  general  form  of  the  response  to  impulsive 
loading  was  evaluated  using  contour  integration. 

In  conclusion,  the  approach  to  viscoelasticity  resulting 
from  this  investigation  is  particularly  powerful,  in  that  the 
general  motion  of  a structure  having  both  elastic  and  visco- 
elastic components  can  be  determined  given  that  two  conditions 
are  met.  First,  the  RT  and/or  RTG  models  for  the  viscoelastic 
materia.ls  in  the  structure  must  exist  and  comply  with  the 
second  law  of  thermodynamics  as  indicated  in  Appendix  A. 

Second,  every  component  of  the  structure  must  have  a suitable 
finite  element  approximation. 


106 


Bibliography 


1.  Caputo,  M.  "Linear  Models  of  Dissipation  Whose  Q Is 
Almost  Frequency  Independent,"  Ann.  Geofisica,  XIX  (4): 
383-393  (1966). 

2.  Christensen,  R.M.  Theory  of  Viscoelasticity,  An 
Introduction.  New  York:  Academic  Press,  1971,  p 14. 

3.  Ross,  B.  "A  Brief  History  and  Exposition  of  the 

Fundamental  Theory  of  Fractional  Calculus,"  Lecture 
Notes  in  Mathematics,  4 57  : 19  . New  York:  Spnnger- 
V'erl'ag,'  19  75  . 

4.  Christensen,  R.M.  Theory  of  Viscoelasticity,  An 
Introduction . New  York:  Academic  Press,  1971,  pp  207-218. 

5.  Hurty,  W.C.  and  M.F.  Rubenstein.  Dynamics  of  Structures, 
Prentice-Hall,  1964,  p 257. 

6.  Crandall,  S.H.  "Dynamic  Response  of  Systems  with 
Structural  Damping,"  Air,  Space  and  Instruments,  Draper 
Anniversary  Volume,  edited  oy  H.S.  Lees . New  York: 
McGraw-Hill,  1963,  pp  183-193. 

7.  Milne,  R.D.  "A  Constructive  Theory  of  Linear  Damping," 
Proceedings  of  Symposium  on  Structural  Dynamics. 
C.4.1-C.4.19.  Loughborough , England : Loughborough 
University  of  Technology,  23-25  March  1970. 

8.  Caputo,  M.  "Vibrations  of  an  Infinite  Plate  with  a 
Frequency  Independent  Q,"  J.  Acoust.  Soc . Am.,  60:  637 
(1976). 


9.  "Linear  Models  of  Dissipation  Whose  Q is  Almost 

Frequency  Independent  - II,"  Geophys . J.  R.  Astr.  Soc, 
1_3 : 529-539  (1967). 

10.  Elasticita  e Dissipazione , Bologna:  Zanichelli, 


1969. 

11.  Caputo,  M.  and  F.  Minardi.  Pure  Appl.  Geophys.,  91: 
134-147  (1971). 

12.  . Riv.  Nuovo  Cimento,  2_ : 161-198  (1971). 

13.  Caputo,  M.  "Vibrations  of  an  Infinite  Viscoelastic 

Layer  with  a Dissipative  Memory,"  J.  Acoust.  Soc.  Am,  56 
(3):  898  (1971).  — 


i 


14. 


T 


I 


) 

► 


i 

1 


"Vibrations  of  an  Infinite  Viscoelastic  Layer 
with  a Dissipative  Memory,"  J.  Acoust.  Soc . Am.,  56  (3): 
897-904  (1971).  " 


15.  "Vibrations  of  an  Infinite  Plate  with  a 

Frequency  Independent  Q,"  J.  Acoust.  Soc.  Am.,  60: 
635-636  (1976). 

16.  . "Vibrations  of  an  Infinite  Plate  with  a 

Frequency  Independent  Q,"  J.  Acoust.  Soc.  Am.,  60: 

634  (1976). 

17.  Foss,  K.A.  "Co-ordinates  Which  Uncouple  the  Equations 
of  Motion  of  Damped  Linear  Dynamic  Systems,"  J.  Appl. 


Mech.  , 2_5_:  361  (1958). 

18.  Caputo,  M.  "Vibrations  of  an  Infinite  Plate  with  a 

Frequency  Independent  Q,"  J.  Acoust.  Soc.  Am.,  60:  902 
(1976).  — 

19.  Ross,  B.  "A  Brief  History  and  Exposition  of  the 
Fundamental  Theory  of  Fractional  Calculus,"  Lecture 
Notes  in  Mathematics,  4 57:  18  . New  York:  Spnnger- 
Verlag,  1975. 

20.  Caputo,  M.  "Vibrations  of  an  Infinite  Plate  with  a 
Frequency  Independent  Q,"  J.  Acoust.  Soc.  Am.,  60:  10 
(1976). 

21.  Gurtin,  M.E.  and  E.  Sternberg.  "On  the  Linear  Theory 
of  Viscoelasticity,"  Arch.  Ration.  Mech.  Anal.  11:  291 
(1962). 

22.  Pipkin,  A.C.  "Approximate  Constitutive  Equations," 

Modern  Developments  in  the  Mechanics  of  Continua, 
Proceedings  of  an  International  Conference  Center  on 
Rheology  held  at  the  Pinebrook  Conference  Center  of 
Syracuse  University,  Pinebrook,  New  York,  23-27  August 
1965.  95.  Edited  by  S.  Eskinazi,  New  York:  Academic 
Press,  1966. 

23.  Christensen,  R.M.  Theory  of  Viscoelasticity , An 
Introduct ion . New  York:  Academic  Press,  19/1  , p 42. 

24.  Zienkiwicz,  O.C.  The  Finite  Element  Method  in  Structural 
and  Continuum  Mechanics.  London:  McGraw-Hill,  1967, 

p 220. 

25.  The  Finite  Element  Method  in  Structural  and 

Continuum  Mechanics.  London:  McGraw-Hill,  1967,  pp  11-12. 


m 


108 


*1 


26.  Foss,  K.A.  "Co-ordinates  Which  Uncouple  the  Equations 
of  Motion  of  Damped  Linear  Dynamic  System^,"  J.  Appl. 

Mech.  , 2_S : 361  (1958)  . 

27.  Bisplinghof  f , R.L.,  et_  al . Aeroelasticity  . Reading, 
Massachusetts:  Add is on -Wes ley , 1955,  p 168. 

28.  Truesdell,  C.  "Thermodynamics  of  Deformations," 

Modern  Developments  in  the  Mechanics  of  Continua, 
Proceedings  of  an  International  Conference  Center  on 
Rheology  held  at  the  Pinebrook  Conference  Center  of 
Syracuse  University,  Pinebrook,  New  York,  23-27  August 
1965.  1.  Edited  by  S.  Eskinazi,  New  York:  Academic 

Press,  1966. 

30.  Coleman,  B.D.  Arch.  Rational  Mech.  Anal.,  17:  1 
(1964  a) 

31.  Donaldson,  J.A.  "A  Family  of  Integral  Representations 
for  the  Solution  of  the  Diffusion  Equation,"  Lecture 
Notes  in  Mathematics,  457:  146.  New  York:  Sprmger- 
Verlag  , l97  5 . 

32.  Coleman,  B.D.  Arch.  Rational  Mech.  Anal.,  1 7 : 230  (1964  b) 

33.  Schlichting,  H.  Boundary  Layer  Theory.  New  York: 
McGraw-Hill,  1968,  p 85. 

34.  Jones,  D.I.G.  "Temperature-Frequency  Dependence  of 
Dynamic  Properties  of  Damping  Materials,"  Journal  of 
Sound  and  Vibration,  33  (4):  451-470  (1974JT 

35.  Christensen,  R.M.  Theory  of  Viscoelasticity,  An 
Introduction . New  York:  Academic  Press,  1971  , pp  96-98. 


l 


4 


109 


Appendix  A 

The  RT  and  RTG  Models  and  the  Second  Law  of  Thermodynamics 

The  focus  of  the  thermodynamic  considerations  centers 
on  the  ability  of  the  RT  and  RTG  models  of  viscoelastic 
materials  to  satisfy  the  second  law  of  thermodynamics.  In 
particular,  the  objective  is  to  derive  constraints  on  the 
parameters  of  the  RT  and  RTG  models  to  insure  that  the  state 
of  stress  and  strain  in  the  material  predicted  by  the  models 
is  consistent  with  the  second  law.  The  considerations  are 
limited  to  the  case  where  the  material  is  at  a constant, 
uniform  temperature  and  the  deformations  and  resulting 
stresses  in  the  material  are  sinusoidal. 

The  assumption  of  constant,  uniform  temperature  is 
motivated  by  the  observation  that  the  moduli  of  viscoelastic 
materials  are  usually  strongly  dependent  on  temperature. 
Allowing  the  temperature  of  the  material  to  vary,  due  to 
the  energy  dissipated  during  the  deformation  process,  implies 
that  the  parameters  of  the  RT  and/or  RTG  model  should  be 
modified  to  account  for  the  resulting  changes  in  material 
properties.  However,  the  RT  and  RTG  models  are  employed 
where  the  parameters  of  the  models  are  not  changed  during 
the  deformation  process,  which  is  in  essence  a statement 
that  the  temperature  of  material  does  not  change  during 
deformation . 


To  implement  the  assumption  of  constant,  uniform 
temperature,  energy  sinks  are  assumed  to  be  present  at 
every  point  in  the  material.  The  energy  sinks  instantane- 
ously absorb  all  the  dissipated  energy,  preventing  the 
temperature  from  changing. 

The  energy  sinks  are  represented  in  the  first  law  of 
thermodynamics 

pe  = p + V*fi  + pq  A.l 

by  q,  the  local  rate  of  change  of  energy  density  through 
non-mechanical  effects.  In  this  local  form  of  the  first 
law,  is  the  divergence  of  the  energy  flux  vector,  p 

is  the  local  rate  at  which  internal  work  is  done  by  mechani- 
cal means,  e is  the  internal  energy  per  unit  mass  of  the 
material,  and  p is  the  mass  density. 

The  following  development,  based  solely  on  the  first 
and  second  laws,  is  an  adaptation  of  the  work  of  Coleman 
and  Truesdell  on  the  thermodynamics  of  deformations  for  the 
specific  case  at  hand  (Ref  29).  In  notation  identical  to 
Truesdell's,  the  second  law  for  a material  at  constant, 
uniform  temperature  is 

pen  - v*h  - pq  > 0 A. 2 

where  n is  the  entropy  per  unit  mass  and  0 is 
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temperature.  This  local  form  of  the  second  law  is  based 
on  the  Clausius -Planck  inequality. 


Two  other  entities  relevant  to  the  following  discussion 
are  the  free  energy  of  the  material,  <|>  , 

<|>  = e - n0  A. 3 

and  the  internal  dissipation  rate,  6 , 

6 = p - p(ne  - i)  A. 4 

where  the  dots  indicate  first  time  derivatives.  Substituting 
the  expression  for  the  dissipation  rate  produces 

6 = pen  + (p  - pe)  A. 5 

Using  the  first  law  to  substitute  for  p - pe  in  Eqn  A. 5 
yields 

6 = pen  * v*h  - pq  A. 6 


In  other  words,  when  the  local  internal  dissipation  rate  of 
a material  is  non-negative,  the  behavior  of  the  material 
is  consistent  with  the  second  law.* 

The  task  remaining  is  to  express  the  dissipation  rate 
in  terms  of  the  stresses  and  strains  in  the  material.  To 
that  end,  recall  that  the  energy  sinks  absorb  all  the 
dissipated  energy. 

6 = - pq  A . 8 

In  light  of  the  expression  for  the  dissipation  rate  given 
in  Eqn  A. 6, 


-v  *FT  + pen  = 0 


A.  9 


Since  the  absorption  of  the  energy  by  the  sinks  prevents 
the  opportunity  for  energy  conduction, 

K = 

and,  as  a result 

7*h  = 0 

* This  result  is  contained  in  Coleman's  general  theorem  on  the  non- 
equilibrium thermodynamics  of  deformation  (Refs  30,  32). 
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and  consequently 


n = 0.  A. 12 

Since  the  rate  of  change  of  the  entropy,  n , is  zero, 

Eqn  A. 5 reduces  to 

6 = (p  - pe)  A. 13 

The  internal  dissipation  rate  can  now  be  expressed  in 
terms  of  stress  and  strain,  provided  that  the  internal 
energy  per  unit  volume,  pe  , can  be  expressed  in  terms 
of  the  rate  of  change  of  strain  energy,  u . The  internal 
energy  per  unit  mass,  e , is  taken  to  be  a function  of 
temperature  plus  the  strain  energy  per  unit  mass,  u^p 

e = f (6)  + u/p  A. 14 

For  small  strain  the  expression  for  one  over  density  is 

17P(  t)  - (1  + e(t))/pQ  A. 15 

where  pQ  is  the  unstrained  density  and  e(t)  is  the 

dilatation  strain,  defined  in  Eqn  4.5.  The  resulting 
expression  for  the  internal  energy  is 
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e = f (6)  + u*  (1  + e (t) )/p 


A. 16 


which  simplifies  to 


e = f (6)  + u/p 


A. 17 


for  small  strain.  Since  the  temperature  is  constant,  the 
expression  for  the  strain  energy  rate  is 


• : u/ 

e p. 


A. 18 


and  the  resulting  expression  for  the  dissipation  rate  is 


4 


6 = (p  - u) 

again,  for  the  small  strain. 

The  dissipation  rate  in  terms  of  the  stresses  and 
strains  is 

6 - ({  o (t  ) }T{  e (t ) } - {o'(t)}T{e(t)}) 


A.  19 


A.  20 


or 


6 = Uo(t)}  - {o'(t)}}l{e(t)} 


where  {e(t)}  is  the  column  vector  of  the  strain  rates, 


A.  21 


— • ~ — ^ — " 
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{ 0 (t ) } is  the  column  vector  of  the  stresses  and  (o'(t)} 
is  the  column  vector  of  those  components  of  the  stress 
acting  to  store  and  retrieve  strain  energy  in  the  material. 
As  indicated  in  Section  IV,  the  stresses  storing  and 
retrieving  strain  energy  for  steady  - state , sinusoidal 
motion  were  those  stresses  in-phase  with  the  strains. 
Consequently,  the  total  stresses,  (o(t)}  , minus  the 

in-phase  stresses,  (o'(t)}  , leaves  the  out-of-phase 

stresses,  {o"(t)} 

6 = {o"(t)}T{e(t)  } 


The  stress  out-of-phase  with  the  strains 


e„„(t)  = e sinu  t 

mn  mn  o 

o 


are  those  stresses  proportional  to  cos^t  in  Eqn  4.32 


a"  ft)  = 6 )e  cosu>  t + F.  (w  )e  „ cosu  t 

mn  mn  3 o'  o o 4vo  mnQ  o 


The  dissipation  rate  in  terms  of  the  stresses  and  strains 
given  above  is 


3 3 

6 = ui  F,  (w  )e2cos2u  t + 7 J u)  F.(u>  )e2  cos2w  t 

o 3V  oJ  o o i o 4 v o mn_  o 


n=l  m-1 


* Fj(a>o)  and  F^(a)Q)  are  functions  of  the  parameters  of  the  RT  and 
RTG  models  and  are  defined  by  Eqns  4.34  and  4.37. 


A. 22 


4.31 


A. 23* 


A. 24 
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If  the  parameters  of  the  RT  and/or  RTG  models  are 
chosen  such  that 

“oF3(“0>  i 0 A. 25 

<1>oF4^wo^  - 0 A. 26 

for  all  positive  uQ  , the  state  of  stress  and  strain  in 
the  material  predicted  by  the  models  is  consistent  with  the 
second  law. 

Although  the  above  development  is  for  sinusoidal  motion 
of  the  material,  its  application  extends  to  periodic  strain 
histories.  Consider  the  case  where  the  material  undergoes 
prescribed,  piece-wise  continuous  and  bounded  strain  history 
of  finite  duration  and  the  material  is  allowed  sufficient 
time  for  the  stresses  to  relax  to  their  unloaded  equilibrium 
value  of  zero.  Assume  that  at  some  later  time  the  same 
strain  history  is  imposed  on  the  material  and  the  material 
again  allowed  to  fully  relax,  and  this  process  of  loading  and 
relaxation  is  repeated  continuously  at  a specified  interval. 

The  stress  and  strain  histories  of  the  material  are  periodic 
and  can  be  broken  into  their  respective  frequency  components 
using  Fourier  analysis.  The  internal  dissipation  of  each 
frequency  component  of  the  strain  is  given  by  Eqn  A. 24  and 
satisfies  the  second  law,  given  that  Eqns  A. 25  and  A. 26  hold. 
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Consequently,  the  periodic  stress  and  strain  histories 
of  the  material  comply  with  the  second  law.  As  a result, 
satisfying  A. 25  and  A. 26  is  sufficient  to  insure  that  the 
response  of  the  material  to  a piece-wise  continuous  and 
bounded  strain  history  of  finite  duration  satisfies  the 
second  law.  This  follows  from  the  observation  that  each 
loading  and  relaxation  cycle  of  the  material  can  be  viewed 
as  an  independent  response  sequence. 


Appendix  B 

A Generalized  Derivative  Relation  for  a Newtonian  Fluid 


By  way  of  motivation  toward  describing  the  motion  of 
physical  systems  with  generalized  derivative  equations,  it 
is  interesting  to  note  that  a Newtonian  fluid,  defined  by 
the  stress -strain  rate  constitutive  relation 


B.  1 


where  p is  the  bulk  viscosity,  can  undergo  dynamic  motion 
which  lends  itself  to  description  using  a generalized 
derivative.  When  the  Newtonian  constitutive  relation  is 
introduced  into  the  one-dimensional  momentum  equation  for 
a homogeneous,  incompressible  fluid  at  uniform,  constant 
temperature 


3 v _ 3 o 

P 3t  = 3 Z 

* 

the  resulting  differential  equation  is 


B.2 


has  demonstrated  that  solutions  of  the  one -dimens ional 
diffusion  equation  may  be  represented  using  fractional 
calculus.  Of  particular  interest  at  this  point  is  the 
solution  to  the  one -dimens ional  diffusion  equation  that 
occurs  when  an  infinite  half-space  of  such  a Newtonian 
fluid  is  bounded  by  a "wetted"  planar  surface  undergoing 
some  known  motion.*  If  the  known  motion  of  the  bounding 
surface  has  a Laplace  transform,  Laplace  transforms  may  be 
used  to  solve  the  diffusion  equation. 

A derivation  of  the  solution  is  given  here  for 
completeness.  Taking  the  Laplace  transform  of  Eqn  B.3 
results  in 

pCsV(s.z)  - v(0,z))  = p (s,z) 

3 Z 2 

where  V(s,z)  is  the  transform  of  the  velocity 

00 

V(s,z)  = / v(t,z)e  stdt 

Assuming  the  fluid  is  initially  excited  from  a state  of 
rest,  the  initial  velocity,  v(0,z)  , is  zero  and  the 

resulting  differential  equation  is 


* This  problem  is  a generalization  of  Stokes'  second  problem  in  which 
Stokes  assumed  the  motion  of  the  "wetted"  plate  to  be  sinusoidal 
(Ref  33). 
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Newtonian  Viscous 
Fluid 


"Wetted”  Surface 


Schematic  of  the  Half  Space  of 
Newtonian  Fluid  Bounded  by  a 
"Wetted"  Surface 


Figure  B. 1 . 


z 


Axis 


psV (s , z ) 


(s  ,z) 


B . 6 


1 


32V 


3 Z ‘ 


The  solution  of  this  differential  equation  is  readily  seen 
to  be 


V (s , z ) 


A(s)e+ 


z 


+ 


B(s)e 


B . 7 


Choosing  the  coordinate  system  such  that  the  normal 
of  the  bounding  surface  is  parallel  to  the  negative  Z axis, 
one  can  proceed  to  determine  the  functions  A(s)  and  B(s) 
in  terms  of  the  boundary  conditions.  The  first  boundary 
condition  to  be  satisfied  is  that  the  magnitude  of  the 
velocity  of  the  fluid  be  bounded  for  large,  negative  Z 
Consequently,  B(s)  must  be  zero.  The  other  boundary 
condition  is  that,  at  the  interface  of  the  fluid  and  the 
bounding  surface,  the  velocity  of  the  fluid  must  be  equal  to 
the  velocity  of  the  bounding  surface,  vQ(t)  • To  enforce 
this  boundary  condition,  one  first  takes  the  transform  of 

vo(t) 


V0(s) 


/ v (t)e'stdt 
0 0 


B.  8 


and  the  boundary  condition  is  applied  in  the  transform  domain 
as  shown 


V (s  , z ) 


A(s)e* 


z 


B . 9 
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V(s  ,0)  = A(s)  = Vo(s) 


B.  10 


The  resulting  expression  for  the  transform  of  the  velocity 
of  the  fluid  at  a distance  z below  the  plate  is 


fps 


V(r,z) 


VQ(s)e 


B.ll 


At  this  point  one  wishes  to  determine  the  stresses  in 
the  fluid  according  to  the  Newtonian  constitutive  relation, 
Eqn  B.l.  The  transform  of  Eqn  B.l  is 


o*(s  ,z) 


9V  (s , z) 
3 Z 


B.  12 


Substituting  Eqn  B.ll  into  Eqn  B.12  for  V(s,z)  produces 


o*(s,z)  = /pp  /s  V(s,z) 


B.  14 


The  transform  of  the  stress  may  be  inverted  by  observing  that 
the  transform  is  the  product  of  the  transforms  of  two  known 
functions  of  time. 


o*(s,z)  = /up  yj  • sV  (s  , z ) = /ppL[  (r  (>j)tli) _1] 


• L[v (t , z ) ] 


B.  15 
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Thus,  the  stress  is  seen  to  be  the  convolution  of  the  two 
functions  of  time, 


1 


< 


a(  t,z)  = dt.  B.  16 

r (‘j)  o (t-xT 

For  the  zero  initial  condition  on  velocity,  Eqn  B.16  is 
equivalent  to 

o(t,z)  = A.  J 1 dt  B.  17 

r(»j)  dt  0 (t-t)* 

Notice  that  Eqn  B.17  states  that  the  stress  at  any 
location  in  the  fluid  is  equal  to  a constant,  /yp  , times 
the  generalized  derivative  of  fractional  order  of  the 
velocity  of  the  fluid  at  that  location. 

o(t,z)  = /y"p  D*5  [v  (t , z ) ] B.  18 

This  stress -velocity  relation  evaluated  over  an  area  A of 
the  "wetted"  surface,  z = 0 , produces  a force -velocity 

relation. 

f (t ,0)  - A JVe  Dh  [v0(t)]  B. 19 


4 


V 
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Thus,  we  see  that  in  this  one  case,  the  macroscopic  behavior 
of  a Newtonian,  viscous  fluid  is  characterized  by  a generalized 
derivative  of  fractional  order  , even  though  the  micro- 

scopic behavior  is 

o(t,z)  = ye(t,z)  B.2 

This  observation  suggests  that  generalized  derivatives  may 
have  applications  in  other  situations  where  a global,  or 
discrete,  description  is  desired  of  a phenomenon  which  is 
locally  viscous. 


Appendix  C 

RT  and  RTG  Models  for  Viscoelastic  Materials 


Generalized  derivative  constitutive  relations  for 
three  viscoelastic  materials  are  presented.  The  material 
properties  on  which  the  RT  and  RTG  models  are  based 
(Ref  34)  were  determined  using  the  "temperature-frequency 
superposition"  principle  (Ref  33). 

The  RTG  model  for  the  shear  modulus  of  the  visco- 
elastic material  3M-467  at  75°F  is 


(1  ♦ bjD  A)o(t)  = (pQ  ♦ pxD  A)e(t) 


where 


bj  = 8 x 10’4  sec*  51 

uQ  = 1.0  lb/in2 

Vj  = 7.3  lb-sec* 56/in2 

- .51 


1 


C.  2 

C.  3 

C . 4 

C . 5 


C.  6 
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The  mechanical  properties  predicted  by  the  model  are 
compared  to  the  material's  properties  in  Fig  C.l.  Note 
the  excellent  agreement  between  the  model  and  the  material 
properties  over  8 decades  of  frequency. 

The  RT  model  for  the  viscoelastic  Young's  modulus  of 
Sylgard  188  at  120°F  is 


o(t)  = (E0  + EXD  )e(t) 


where 


Eq  = 60  lb/in2 


24 


Ej  = 43  lb-sec*  /in2 


and 


.24 


The  mechanical  properties  predicted  by  the  model  are 
compared  to  the  material's  properties  in  Fig  C.2. 

The  RT  model  for  the  viscoelastic  Young's  modulus 
of  BTR  at  45°F  is 


o(t)  - (E^"1  ♦ E2D°2)e(t) 


C . 7 


C.  8 


C . 9 


C.10 


C.ll 
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where 


C.  1 2 

C.  13 

= . 095  C . 14 

and 

c*2  = *28  C . 1 5 

A comparison  of  the  model  and  material  properties  is  given 
in  Fig  C . 3 . 

Although  the  agreement  between  the  material  properties 
and  their  respective  models  is  very  good,  not  all  visco- 
elastic materials  lend  themselves  to  characterization  by 
generalized  derivatives  models.  The  materials  most  suited 
to  modeling  with  generalized  derivative  constitutive 
relations  are  those  that  have  properties  such  that  the 
following  relation  holds  in  the  transition  region: 

n = tan  — C.16 

2 

where  n is  the  loss  factor  and  a is  the  slope  of  the  plot 
of  the  log^o  of  the  real  part  of  the  modulus  plotted  as  a 
function  of  log1Q  of  the  frequency  of  motion. 


09  5 

= 8S0  lb-sec’  /in2 

E-  = 18  lb-sec • 28/in2 
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In  summary,  generalized  derivative  constitutive  relations 
do  in  fact  model  the  frequency  dependent  mechanical  properties 
of  at  least  three  viscoelastic  materials.  However,  each 
model  is  of  a different  basic  form,  as  can  be  seen  by  comparing 
Eqns  C.l,  C.7,  and  C.ll.  Hence,  the  RT  or  RTG  models  are 
capable  of  describing  viscoelastic  materials  having  distinctly 
different  mechanical  properties. 


The  Mechanical  Properties  of  3M-467 
Compared  to  the  RTG  Model  of  3M-467 


Figure  C.2.  The  Mechanical  Properties  of  Sylgard  188 
Compared  to  the  RT  Model  for  Sylgard  188 
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Generalized  derivatives  of  fractional  order  are  used  to  construct  stress- 
strain  constitutive  relations  for  viscoelastic  materials,  based  on  the 
observed  sinusoidal  behavior  of  the  materials.  The  non-periodic  behavior 
of  one  material  is  observed  in  the  laboratory  and  compares  favorably  with 
the  non-periodic  behavior  of  the  material  predicted  by  its  generalized 
derivative  constitutive  relation.  Having  established  that  the  generalized 
derivative  constitutive  relation  is  an  appropriate  mathematical  model  for 
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the  general  motion  of  at  least  one  viscoelastic  material,  the  tools  for 
the  analysis  of  structures  of  engineering  interest  are  put  forward.  In 
particular,  attention  is  focused  on  a finite  element  formulation  of  and 
solutions  to  the  equations  of  motion  for  structures  containing  elastic 
and  viscoelastic  components. 
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